First-principles model potentials for lattice-dynamical studies: general methodology 
and example of application to ferroic perovskite oxides 
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We present a scheme to construct model potentials, with parameters computed from first princi- 
ples, for large-scale lattice-dynamical simulations of materials. Our method mimics the traditional 
solid-state approach to the investigation of vibrational spectra, i.e., we start from a suitably chosen 
reference configuration of the material and describe its energy as a function of arbitrary atomic 
distortions by means of a Taylor series. Such a form of the potential-energy surface is completely 
general, trivial to formulate for any compound, and physically transparent. Further, the approxi- 
mations involved in our effective models - i.e., the truncations affecting the order of the polynomial 
expansion, the spatial range of the interatomic couplings, and the maximum number of atoms (or 
bodies) involved in the interaction terms of the series - are clear-cut, and the precision can be im- 
proved in a systematic and well-defined fashion. Moreover, such a simple definition allows for a 
straightforward determination of the parameters in the low-order terms of the series, as they are 
the direct result of density-functional-perturbation-theory calculations, which greatly simplifies the 
model construction. Here we present such a scheme, discuss a practical and versatile methodology 
for the calculation of the model parameters from first principles, and describe our results for two 
challenging cases in which the model potential is strongly anharmonic, namely, ferroic perovskite 
oxides PbTiOs and SrTi03. The choice of test materials was partly motivated by historical reasons, 
since our scheme can be viewed as a natural extension of (and was initially inspired by) the so-called 
first-principles effective Hamiltonian approach to the investigation of temperature-driven effects in 
ferroelectric perovskite oxides. Thus, the study of these compounds allows us to better describe the 
connections between the effective-Hamiltonian method and ours. 

PACS numbers: 63.70.+h,71.15.Mb,65.40.-b 



I. INTRODUCTION 

The development of methods for statistical simula- 
tions with first-principles accuracy remains one of the 
major challenges for the community working on com- 
putational condensed-matter physics and materials sci- 
ence. In spite of recent advances, state-of-the-art first- 
principles methods are still unable to reach the length 
and time scales that are relevant for the study of many 
properties of interest at realistic operating conditions. 
Ranging from temperature-driven phase transitions to 
thermally-activated processes of all sorts, there are count- 
less phenomena whose first-principles treatment has a 
prohibitive computational cost, even if one resorts to 
the most numerically-efficient schemes such as density- 
functional theory (DFT)ji Hence, there is a need to de- 
velop approximate methods that allow for fast calcula- 
tions while retaining the first-principles accuracy and, if 
possible, predictive power. Much of the on-going activ- 
ity on multi-scale simulations is the direct consequence 
of this situation. 

Whenever one is concerned with the lattice-dynamical 
properties of the materials, it may be possible to avoid 
the explicit treatment of the electrons in the simulations. 
Such is typically the case when we are interested in struc- 
tural and mechanical properties, dielectric and piezoelec- 
tric responses (which are dominated by the lattice part of 
the effect, as opposed to the electronic one, in the mate- 



rials that are most attractive for applications), or lattice 
thermal transport, to name a few important examples. 
Many methods have been developed to address this sub- 
set of problems, which are the focus of the present work. 

In the context of lattice-dynamical studies, there are 
essentially two families of effective potentials that allow 
for large-scale simulations. The most widely used meth- 
ods have been developed in the Physics and Chemistry 
communities. Such models - which include Lennard- 
Jones potentials^ shell models^ bond- valence models^ 
and even Tersoff potentials^ and reactive force fields? 7 - to 
name a few - tend to rely on a physically-motivated ana- 
lytic form of the atomic interactions. Unfortunately, that 
restriction is often a too stringent one, and compromises 
the ability of the models to reproduce the first-principles 
data. Further, a systematic extension of the models to 
improve precision is usually not well defined or possible. 

Another approach is represented by the methods that 
rely on artificial neural networks^ importing techniques 
developed by the artificial intelligence community. In 
this case, one uses very versatile models that can repro- 
duce first-principles data with arbitrary precision, at the 
expense of creating complicated potentials that do not al- 
low for a clear physical interpretation. On top of the loss 
in fundamental understanding, the fact that such models 
are not physically motivated usually implies that they 
have poor transferability and a limited predictive power 
[i.e., they are good for interpolating between the first- 
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principles data points used to fit the model, but often 
fail when used to predict (extrapolate) new behaviors]. 
Additionally, they are relatively costly from a computa- 
tional point of view, as the potentials can become quite 
complex. 

Interestingly, some workers have developed alterna- 
tive, very successful approaches that overcome most of 
the above mentioned deficiencies, but which have been 
applied to a very small set of problems. One relevant 
example is the work of Rabe, Vanderbilt, and others 
on ferroelectric perovskite oxides: Already in the 1990s 
these authors constructed first-principles model poten- 
tials, which arc usually called effective Hamiltonians, 
to investigate the ferroelectric phase transitions of per- 
ovskites like BaTiO^^ and PbTi0 3 .-±i [These works 
built upon the ideas introduced in Ref. Qj] to investi- 
gate the structural phase transition in GeTe from first 
principles.] The effective-Hamiltonian approach involves 
a drastic simplification of the material, which is coarse- 
grained to retain only those degrees of freedom associ- 
ated with the ferroelectric properties (i.e., local dipolcs 
and cell strains). The potential-energy surface (PES) 
corresponding to these relevant variables is written as 
a low-order Taylor series around a suitably chosen refer- 
ence structure (i.e., the prototype cubic perovskite struc- 
ture). Such a scheme is physically- motivated, compu- 
tationally very efficient, and its precision can be im- 
proved, to some extent, in a well-defined way. Further, 
the application of the original approach to increasingly 
complex oxides (e.g., compounds with non-polar tran- 
sitions like SrTi03 ) 13 i 14 chemically-disordered materials 
like PbZri-^TirOa ) 15 ^ 6 and magnetoelectric multifer- 
roics like BiFe0 ^ 17 i 18 ) has shown its generality, the good 
transferability of the interatomic couplings among dis- 
similar chemical environments, and the reliability and 
predictive power of the models. Unfortunately, as far as 
we know, such an approach has not been adopted in other 
research fields, remaining much confined within a small 
community working on ferroic perovskite oxides. 

In our opinion, to understand why the effective- 
Hamiltonian method has failed to gain widespread pop- 
ularity, one has to consider the coarse-graining step in- 
volved in the construction of the potential. When these 
models were first developed, there were plenty of reasons 
to adopt such a simplification. On the one hand, by re- 
stricting to a subset of the configuration space, it is pos- 
sible to construct simpler potentials and run faster sim- 
ulations. On the other hand, by the time first-principles 
methods started to be applied to these problems, there 
was already a whole body of literature devoted to simi- 
lar, semi-empirical models used in theoretical studies of 
phase transitions driven by soft modes* 1 ^— Indeed, the 
effective Hamiltonians of Refs. and [111 can be viewed 
as the natural evolution of the models that already ex- 
isted in the literature, as for example the so-called dis- 
crete <j) A models In some sense, the main innovation in 
those pioneering works was to develop a systematic and 
well-defined scheme to compute the parameters of such 



Hamiltonians from first principles. To do that, the key 
step was to establish a connection between the variables 
of the traditional effective models (i.e., the so-called lo- 
cal modes representing the localized atomic distortions 
whose collective occurrence leads to the structural tran- 
sition, and which involve the formation of local electric 
dipoles in the case of fcrroelectrics) and the displace- 
ments of the actual atoms in the crystal. Such a connec- 
tion can be made in a variety of ways, ranging from the 
more elementary 1 ^ (e.g., by defining the local modes from 
direct inspection of the strongest structural instabilities 
of the high-symmetry phase, which can be determined 
from first principles as discussed below) to the more 
sophisticate d 23 ^ 4 (e.g., by identifying the local modes 
with lattice Wannicr functions computed from knowledge 
of the full phonon dispersion bands of the high-symmetry 
phase as obtained from first principles). Once the local 
modes are defined in terms of actual atomic displace- 
ments, the first-principles calculation of the Hamiltonian 
parameters follows in a rather straightforward way. 

However, in general there are no clear reasons to in- 
troduce such a coarse graining. Suppose, for example, 
that we want to investigate a stable phase of a mate- 
rial, and need a model that captures the first-principles 
energetics with very high precision. In such a case, in 
absence of structural instabilities of our reference con- 
figuration, it may be unclear how to choose a subset 
of relevant degrees of freedom. Further, we may typ- 
ically find that most of the modes, even the relatively 
high-energy ones, play an important role in determin- 
ing the properties of such a phase at a quantitative level; 
hence, the restriction to a configuration subspace, and 
the reduced accuracy it entails, may not be acceptable. 
Even in cases where the focus is on the investigation of 
structural phase transitions, and assuming that we would 
be satisfied by a sound qualitative description of such a 
drastic effect, the coarse-graining step may turn out to 
be both unhelpful and difficult to implement. Consider, 
for example, the modern perovskite oxide super-lattices 
that present a wealth of appealing physical effects, some 
of which are attributed to novel interatomic couplings 
occurring at the interfaces between different layers. In 
such situations, suitably exemplified by the short-period 
PbTiOa/Si'TiOs super-lattices studied by some of us^ 
identifying a small set of relevant degrees of freedom may 
be very hard; indeed, a large number of local modes may 
need to be considered, which would result in complicated 
effective models and a relatively small gain in computa- 
tional efficiency. Other cases where similar difficulties 
are likely to appear are crystals in which the relative sta- 
bility of different phases depends strongly on secondary 
structural order parameters (as in the case of BiFcOg 2 ^) , 
situations in which strong strain gradients and non-trivial 
structural relaxations occur (as in the vicinity of ferro- 
electric domain walls 2 ^), etc. In our opinion, the local- 
mode approximation is not well suited for the study of 
such problems. 

In view of this, we decided to adopt an approach 



3 



that retains many of the good features of the effective- 
Hamiltonian method developed within the ferroelectrics 
community and avoids its most serious limitations. In 
short, we decided to create models that describe the en- 
ergetics of all the atomic degrees of freedom of a ma- 
terial by Taylor expanding the PES around a suitably 
chosen reference structure. Working with a simple poly- 
nomial model has many important advantages: (i) it is 
general and can be trivially formulated for any material; 
(ii) the involved parameters have an obvious and con- 
venient physical interpretation, as we essentially follow 
the approach adopted in solid-state textbooks to discuss 
lattice-dynamical and elastic properties; (Hi) many of the 
potential parameters can be obtained directly from per- 
turbative first-principles calculations; and (iv) the pre- 
cision of the potential can be improved in a systematic 
and well-defined way. Here we describe the details of 
such a scheme and illustrate it with applications to two 
challenging cases: ferroics PbTiC>3 and SrTi03, both of 
which present soft-mode-driven structural phase transi- 
tions whose description requires the use of strongly an- 
harmonic potentials. We are thus introducing a method 
that we think should be very useful and of general appli- 
cability. 

The paper is organized as follows. In Section II we de- 
scribe the general methodology, using perovskite oxides 
as an illustrative sample case. We also introduce the ap- 
proach we adopted to compute the potential parameters 
from first principles. In Section III we describe the mod- 
els constructed for perovskite oxides PbTiC>3 and SrTiC>3. 
We also solve the models by means of Monte Carlo simu- 
lations, showing that they reproduce the experimentally 
observed phase transitions. Finally, in Section IV we 
summarize and conclude the paper. 



II. MODEL CONSTRUCTION 

In the following we present our general scheme for con- 
structing effective model potentials for lattice-dynamical 
studies. The proposed methodology is general and can 
be applied to any material, including cases of reduced 
dimensionality (e.g., surfaces, slabs, wires, molecules), 
disordered systems, etc. Nevertheless, for the sake of 
clarity, here we will refer to the case of an infinite peri- 
odic crystal, and take the family of perovskite oxides as 
a representative example of application. 



A. Reference structure and model variables 

The construction of our models begins with the choice 
of a reference structure (RS) that will typically be a min- 
imum or a saddle point of the PES. Thus, for example, 
if we were interested in the properties of a particular 
(meta)stable phase of a material, the RS would corre- 
spond to the solution obtained for such a phase by per- 
forming a first-principles structural relaxation nominally 



at T = K. If we were interested in the more challenging 
case of a material undergoing structural phase transitions 
driven by the softening (i.e., destabilization) of some vi- 
brational modes or cell strains, it would be convenient to 
take the high-symmetry phase of the material as our RS. 
More specifically, in that case we would determine the 
RS by performing a constrained relaxation (i.e., one in 
which the high symmetry of the undistorted phase is im- 
posed) using first-principles calculations at T = K; the 
result will typically be a saddle point of the PES, and will 
have an associated Hessian matrix (i.e., a matrix of sec- 
ond derivatives of the energy) with negative eigenvalues 
corresponding to the structural instabilities. Our cho- 
sen examples of application - i.e., ferroics PbTiC>3 and 
SrTiC>3 - belong to this second category. 

For the sake of concreteness, let us assume that we are 
treating a three-dimensional infinite crystal composed of 
periodically repeated cells. Then, our RS is fully speci- 
fied by the lattice vectors Ri, where I labels cells, and the 
positions t k of the atoms k inside the cell. We will de- 
scribe all accessible crystal configurations as distortions 
of the RS. The most general atomic state will be given 
by the position vectors 

X) (<W + Vaf})(Rlf) + T K p) + U lKa , (1) 



where a and j3 denote Cartesian directions. The distor- 
tions are thus captured by the homogeneous strain r\ and 
the individual atomic displacements u\ K . It is important 
to note that the homogeneous strain affects both the Ri 
and t k vectors defining the RS, and that we describe the 
deviation from the strained RS by means of the absolute 
displacements it; K . Hence, the it; K vectors are given in 
Cartesian coordinates and have units of length. Alterna- 
tively, we could have worked with atomic displacements 
given in units relative to the (strained) reference struc- 
ture; however, our definition of the ui K variables leads 
to a clear and computationally convenient formulation of 
the model potentials, which is why we adopted it. 

The homogeneous strain r] a p in Eq. |T]) contains both 
symmetric and anti-symmetric parts. Typically, it will be 
convenient to exclude the anti-symmetric part (i.e., rigid 
rotations of the whole material) from the description. We 
will thus restrict ourselves to the symmetric components 
(?7q^ +?7^ a )/2, for which we will use the well-known Voigt 
notation ?y a , with a = 1, 6^ 

To alleviate the notation in the formulas below, it is 
convenient to introduce the following bijective mapping 

koi, (2) 

so that any Ik pair can be expressed by a single index i, 
and vice versa. Hence, we can use iti instead of ui K , and 
even write Ri or r.;, without any ambiguity. 

Figure [T] shows the cubic phase of an ABO3 perovskite 
oxide, which is the RS of the applications discussed be- 
low. The shown cell is repeated along the three spatial 
directions and, while the displacements it; may change 
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FIG. 1. Definition of the cubic reference structure for ABO3 
perovskite oxides. The unit cell vectors are a = ao(l,0, 0), b 
= ao(0, 1, 0), and c = ao(0, 0, 1), as expressed in the Cartesian 
reference depicted in the figure. Lattice vectors are given by 
Ri — lino + niib + Tine, where nn, nn, and nn are integers. 
The positions of the atoms within the unit cell are: ta = 
a (0,0,0), t b = a (l/2,l/2, 1/2), roi = a (0, 1/2, 1/2), r 02 
= a (l/2, 0, 1/2), and r 03 = a (l/2, 1/2, 0). 

from cell to cell, the homogeneous strain 77 is the same 
throughout the material. 

In Fig. [5] we illustrate how an arbitrary distortion is 
captured by the variables defined in Eq. ([T]). Panel (a) 
shows an unstrained configuration (r/ = 0) with atoms 
fluctuating around their RS positions. Panel (b) shows 
an homogeneously strained state with the atoms main- 
taining their relative positions {ui = 0). Finally, pan- 
els (c) and (d) show configurations in which, while the 
homogeneous strain is zero, we do have local strains re- 
sulting in the expansion along the vertical direction of 
some cells of the material [(c)] or a shear-like deformation 
[(d)] . Note that such inhomogeneous strains are naturally 
described by the Ui variables in our model. 

Before we continue, it is worth mentioning the advan- 
tages and limitations of constructing model potentials 
based on a reference structure. As will be obvious below, 
this approach makes it trivial to formulate the poten- 
tial for arbitrary materials and with arbitrary precision. 
Hence, it allows for a general and clear-cut methodol- 
ogy that offers the possibility of improving the models 
systematically; those are obviously very important as- 
sets that are not so frequently found in model-potential 
schemes. At the same time, our approach is specifically 
suited for the description of relatively small distortions 
of the RS. In other words, we expect our effective mod- 
els to describe the energy of configurations that resemble 
the RS in some fundamental way; for example, the bond 
topology should be roughly preserved, and we would not 
advice the use of the present method to describe situa- 



tions in which chemical bonds break or form. 

Finally, let us stress that our scheme is applicable to 
any material, not only to the infinite three-dimensional 
crystals that we focus on for the sake of the presenta- 
tion. Indeed, materials of arbitrary dimensionality, or 
disordered compounds, can be tackled by making the ap- 
propriate adjustments. For example, to study a molecule 
we would work with atomic positions defined by r Ka = 
t ko + u Ka ; to work with materials that are periodic only 
along one or two directions, we would just need to con- 
sider an appropriate homogeneous strain tensor; to deal 
with chemically-disordered materials, we would construct 
potentials that depend on the chemical environment of 
the atoms, etc. While some situations may be more chal- 
lenging than others, our methodology can in principle be 
applied in all cases. 

B. Definition of the effective potential 

Based on the variables defined previously, we write the 
energy changes around the RS, AE c g = E c h — £rs, as 

AE eS ({ Ui }, v ) = E p ({ Ui }) + E S ( V ) + E sp ({ Ui }, n), (3) 

where 

E p ({ Ui }) = E ha , T ({ Ui }) + E anh ({ Ul }) . (4) 

Here we use the subscript "eff" to distinguish between 
the energy that our effective potential gives for config- 
uration {{ui],r}) and the real energy E({ui},rj) that 
we would obtain from a first-principles simulation of the 
same atomic state. The above terms are: the energy of 
the RS (-Ers); the energy change when the RS is dis- 
torted by atomic displacements {uf\ (E p , where the "p" 
subscript stands for "phonon"), which we split into har- 
monic (Ehar) and anharmonic (-E an h) contributions; the 
energy change when we strain the RS (E s , where the 
"s" subscript stands for "strain"); and the additional en- 
ergy variations occurring when homogeneous strains and 
atomic distortions appear simultaneously (E sp , where 
"sp" stands for "strain-phonon" ) . Let us discuss each 
of these terms. 

1- -Ehar({ui}) and i? an h({wi}) 

Traditionally, the energy change caused by atomic dis- 
tortions {iii} is written as a Taylor series around the RS 
in the following way 

Ep({Ui}) = 2 H K ll]/3 U ^cUjfS + 

1 ^ (3) (5) 

where the first line shows the harmonic terms included 
in i?har({wi}) and the second line gathers all the higher 
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FIG. 2. (Color on line.) Sketch of various types of distortions and how the associated energy changes are captured by our 
model potentials. Panel (a): atomic displacements in absence of homogeneous strain. Panel (b): homogeneous strain in absence 
of atomic displacements. Panels (c) and (d): local (inhomogeneous) strain rj loc given by atomic displacement patterns; note 
that the distortions in the unstrained areas are rigid translations, which do not contribute to the energy because of the ASR 
satisfied by E p ({ui}) (see the text). 



order terms contained in E an h({v,i}). The tensor 
is formed by the n-th derivatives of the energy, with 



K 



(n) 



d n E eS 



du ia du jl3 . 



IIS 



(6) 



Note that we assume that the RS is a stationary point of 
the PES (i.e., a minimum or a saddle), so that K {1) = 0. 

It is important to realize that the coefficients 
in Eq. ([3]) are not independent. At each order in the 
Taylor series, they are related by the point and lattice- 
translational symmetries of the RS structure. Addition- 
ally, and more fundamentally, they have to comply with 
translational invariance in free space, which results in 
the so-called acoustic sum rules (ASRs). In essence, the 
ASRs guarantee that a rigid translation of the material 
- i.e., one given by Ui a = u a , where u a is an arbitrary 
three-dimensional vector - does not change the energy 
and does not induce any forces on the atoms. To fulfill 
these conditions, the coefficients must satisfy 



/ ^ iajphy... 



0, Vj,fc, 7, (7) 



at all orders of the expansion. In the harmonic case with 
n = 2, this reduces to the well-known ASR for the ele- 
ments of the force-constant matrix 



(8) 



This set of conditions for the harmonic terms is rather 
manageable, and allows for simple procedures to enforce 



the ASR in practice. For example, a common strategy is 
to derive the self-energy parameters from the interactions 
between different atoms, by taking 



K ^ -_V/T (2) 



(9) 



and simultaneously imposing the symmetric character of 
the force-constant matrix {K\J^ = Kjp ia )- Note that 
such a correction is necessary whenever we spatially trun- 
cate the interatomic couplings, as such an approximation 
will generally break the ASR. Also, it is customary to use 
this type of correction when dealing with a force-constant 
matrix whose coefficients may suffer from some numeri- 
cal noise or inaccuracy. As we will discuss in Section lll CI 
we count with well-established and widely-available first- 
principles methods to compute a force-constant matrix 
that is ASR-compliant. Hence, we use the above form 
[i.e., Eq. ([5])] for the harmonic term -Ehar in our models. 

However, as one can imagine from Eq. ([7]), enforcing 
the ASR becomes much more intricate for n > 2. In par- 
ticular, it would complicate enormously the procedure to 
compute the parameters in i? an h discussed in Scction fll CI 
Fortunately, in that case we can resort to an alternative 
representation in which the ASR is automatically satis- 
fied at all orders. 

Indeed, the energy E p ({ui}) can be equivalently ex- 
panded as a function of displacement differences in the 
following way 
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E p({ U i}) =\ K ijLhp( U «* - u jc)( u kP -U h p) + 



ijkh 
a/3 



l E K ^lkh0rt 1 ( U io<-U ja ){u k p-U h p)(u r ~ f -U tl ) + .. 



(10) 



Q / j ' ijotkh/3 
ijkhrt 

a/37 



r 



From this expression, it is obvious that E p does not 
change for a rigid displacement of the material, as ev- 
ery single term cancels out in that case; it is also easy 
to prove that a rigid displacement does not induce any 

~(n) 

forces on the atoms. Hence, the model parameters K 
do not need to satisfy any ASR to guarantee transla- 
tional invariancc, which facilitates enormously the task 
of fitting their values to best reproduce first-principles 
results^ 

The relation between Eqs. (J5J) and (|TU)) is a subtle one 
and deserves a few comments, (i) It is important to 
realize that these two expressions for E p are not con- 
nected by a simple transformation of the basis in which 
we express the atomic distortions of the material. Indeed, 
the atomic displacements {ui} do define the independent 
variables of our problem. In contrast, the set of differ- 



ences {(i 



j a )} has many more, linearly-dependent 



members; hence, the displacement differences are not an 
acceptable basis, (ii) It is possible to go from Eq. ([5]) to 
Eq. (fTU)) by application of the ASR at each order of the 
expansion. More precisely, at a given order n, one can 
use the corresponding ASRs to write some of the K^-^ 
parameters as a function of the rest, e.g. by perform- 
ing substitutions as the one given in Eq. © for n = 2. 
The result of such a procedure is an expression in terms 
of differences, as the one in Eq. (fl0| . However, there 
is no unique way to perform such a transformation and, 
thus, the form of the resulting energy function is some- 
what arbitrary. Indeed, there are many ways in which 
we can use the ASRs to rewrite Eq. ([5]) [e.g., for n = 2, 
Eq. ^ is just one possibility among many others], and we 
see no clear reasons to prefer any specific strategy, (iii) 
In Eq. (|10p. it may look like we have many-body terms 
already at very low orders of the expansion. For exam- 
ple, the harmonic terms can involve up to four different 
atoms. Such couplings are the result of the ASR-related 
connections between the energy derivatives K^ n \ which 
result in terms that look like many-body ones when we 



write the energy as a function of displacement differences. 
(iv) It is possible to understand better the inner struc- 
ture of the difference terms of Eq. (TTU|) . Thus, for ex- 
ample, for n = 2 it can be seen that all four-body terms 
can be written as combinations of two- and three-body 
terms, but three-body terms are in general not reducible 
to two-body terms. These considerations are of little im- 
portance for our present purposes, though, and we will 
not pursue them further, (v) Finally, let us note that the 
~ (") 

K parameters in Eq. (fTOj) can be viewed some sort of 
generalized spring constants; this interpretation is espe- 
cially apparent for the pairwise terms involving products 
of the form (v,i a — Uj a ) n . 



2. E s (r)) and E sp ({ui},T]) 

For the elastic energy E s (r]), we use a simple Taylor 
series 

E M 4S C ab^b + f E C ^aV»Vc + - , (11) 



where 



(m) 



1 d m E eS 



N dr) a dr) b ... 



RS 



(12) 



and TV is the number of cells in the crystal. There is no 
linear term in Eq. (fTTj) because we assume that the RS 
is a stationary point of the PES. The harmonic param- 
eters in this series are the usual elastic constants; more 
precisely, they are the so-called frozen-ion or undressed 
elastic constants, as they quantify the elastic response of 
the material with the ions clamped at the relative posi- 
tions that they have in the RS. 

For the strain-phonon interaction energy E sp ({ui}, tj), 
we can write 



2 Z^i 2-~i aia 'l aUl 



(1,2) 



HP 



1 



V^ A (2,1) 



abia 1 l a7 ^ bUia 



(13) 



The lowest-order coupling term A*- 1 ' 1 -* corresponds (ex- 
cept for non-essential prcfactors) to the so-called force- 
response internal strain tensor, and describes the forces 



I 

that act on the atoms as a consequence of homogeneous 
strains. Hence, this kind of coupling contributes to deter- 
mine the full, relaxed-ion or dressed, elastic response of 
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the material, in the way that is described e.g. in Ref. HO- 
The A < ' m '™' ) parameters in Eq. (p~3|) have to comply with 



with A^ ' parameters that are free from ASR-related 
restrictions. Our choosing between the former or the 
latter expressions for E sp will be a matter of practical 
convenience; more precisely, we will use the regular rep- 
resentation [Eq. (|13p ] whenever we compute the parame- 
ters directly from first principles, and the alternative one 
[Eq. (fT4)) ] in cases in which we need to fit the parameters 
to reproduce specific first-principles results. This will be 
discussed in detail in Section Til CI 



3. Symmetry considerations 

We will often deal with reference structures that 
present certain latticc-translational and/or point symme- 
tries. Such symmetries imply a reduction in the number 
of independent parameters of the model, and we can take 
advantage of them to simplify its construction. In the 
following we describe the general ideas and procedures 
that one can use to this end, resorting to the ABO3 per- 
ovskitcs as a convenient example in which the high sym- 
metry of the RS (i.e., the full cubic space group Pm3m) 
results in great simplifications. 

Let us denote a general symmetry operation by {S^i}, 
where S is the 3x 3-matrix representation of a point sym- 
metry and t is a three-dimensional vector, both expressed 
in our Cartesian reference. By applying such an opera- 
tion to an arbitrary vector x, we obtain the transformed 
vector x' = {S\t}x given by 

X' a = ^2 S al 3Xp + t a . (15) 

p 

For {S^i} to be a symmetry of the RS, it is necessary 
and sufficient to have 

{S\t}(R i + T i ) = R i ,+T i ,, (16) 

where, for any atom i, there is an atom i' of the same 
atomic species that satisfies this relation. In other words, 
Eq. (p~6|) states that the RS is invariant upon the appli- 
cation of {S^i}. 



a set of ASRs that are analogous to the ones discussed 
above for the K^ n ' coefficients. As in the case of E p , we 
can use an alternative expression for -E S p, namely 



(14) 



I 

The distortions of the RS transform as 

U 'i'a=^2 S al3Uip , (17) 

@ 

where i and i' arc related by Eq. (|16[) . and we also have 
vLp = X! Sa~{'n~<s(S~ 1 )8p = ^2 S ai SpsriiS , (18) 

where the strains are expressed in the Cartesian basis. 
Finally, the symmetry condition for the energy reads 

A£off(W},»/') = AE eS ({ Ui },rj) . (19) 

Of course, similar relations hold for all the individual 
terms in the energy and at all orders of the Taylor series 
[e.g., we have E anh ({«•/}) = E anh ({ui})]. 

Let us describe how these general symmetry relations 
allow us to simplify our model potential. Given a par- 
ticular product of Ui a displacements and r\ a strains in 
the Taylor series, we can use the operations of the space 
group to generate the collection of symmetry-related 
products, which will involve transformed u' ila and T]' a 
distortions. Figure [3] illustrates this process pictori- 
ally. For example, in the left panel we start with the 
product {uqbx - uqozx) 2 (u 0By - uoo3y) 2 that couples 
atoms B and 03 in the unit cell at the origin of our 
coordinate system [i.e., without loss of generality, we 
choose Ri = Ro = (0,0,0)]. Then, by application of 
the symmetry operations of the cubic space group of 
the ideal pcrovskite structure, we can generate a col- 
lection of related products; for example, a 90° rotation 
about the y axis transforms the original product into 
(u Bz - u 0O i z ) 2 (u By - uooiy) 2 , etc. Figure [3] sketches 
the products thus generated and involving the B atom 
in the I = cell; lattice-translational symmetry leads to 
analogous couplings centered at all other B atoms in the 
crystal. 

Naturally, these symmetry-related couplings must con- 
tribute to the energy in a very specific way. Continuing 
with the above example, the couplings represented by 
{uqbx - u 0O 3x) 2 (u By - u o3y) 2 appear in our potential 
in the form 



E sp ({Ui}, V) = ^^2Y1 A< aillVa{Uia ~ Uj a ) + J ^ ^flkhpVa^ ~ U Ja )(u k0 - U h[j ) 



a i j a 



a ijhk 
a0 



l Yl ^aUjaVaVbiUia, - U ja ) + 



ab ijoi 
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FIG. 3. (Color on line.) Example of symmetry adapted terms (SATs) involving displacements of neighboring B and 
atoms. We show the SATs represented by the (B x — 03 x ) 2 (B y —03 y ) 2 (left) and (B z — Ol z ) 2 (B x — 01 x ) (right) coupling terms 
(the representative terms are colored in the figure). For the atomic displacements we use the compact notation described 
in the text. An arrow along the a direction, and located at the center of the line connecting atoms i and j, represents the 
(uia — Uj a ) displacement difference. Whenever a displacement difference appears squared, we draw a double arrow, which 
indicates invariance under mirror-plane reflections. Without loss of generality, we assume that the B atom is located at the 
I — cell of the lattice [i.e., Ri — Ro = (0, 0, 0)]; superscripts at the oxygen sites denote the cell to which they correspond. 



K {4) 

"0B0O3;r,0B0O3:i;,0B0O32/,0B0O32/ 



K #15 (B x -O3» ) a (B I( 



(B a 



i u 0Bx — u 0O3x) 2 ( u 0Bi/ — u 0O3y) 2 + ( u 0Bz — u 0Olz) 2 { u 0By ~ u 0Oly) 2 + ■■ 

03? 00 ) 2 + (B a - O3° x 01 ) 2 (B y - O3° 01 ) 2 + (Ba - O2° 00 ) 2 (B, - O2° 00 ) 2 
(B„ ~ 01™ 



(20) 



O2^ 10 ) 2 (B z - 02° 10 ) 2 + (B, - Ol° 00 ) 2 (B y - Ol" 00 ) 2 + (B 2 - Olf °) 2 (B, - OlJ 00 ) 2 



where if #15 is the name we use for this specific coupling 
parameter in Tabled] Here we have introduced a compact 
notation, so that we denote by B a the displacements mobo 
of the central B atom, by Ol™' 1 ™' 2 ™' 3 the displacements 
uioia of the 01 atom at cell /, etc. (Ri is defined by the 
integers nn , 1112, and 1113 as in the caption of Fig. [I]) This 
is what we call a symmetry- adapted term (SAT), which 
is fully specified by one representative coupling and its 
associated parameter. 

Note also that, in some cases, by applying all the space 
group operations to a representative coupling we may 
generate a SAT that exactly cancels out. Hence, work- 
ing with SATs provides us with an automatic way to 
identify couplings that are forbidden by symmetry, which 
may result in drastic simplifications of our model poten- 
tials. For example, in the case of our ABO3 perovskitcs, 
the symmetry of the RS guarantees that all the bilinear 

strain-phonon couplings are zero (i.e., A^ 1 ' 1 ^ = A*' ' = 
0), a fact that has actual physical consequences on the 
response properties of the cubic phase of such materials. 

In our investigation of ABO3 perovskites, we always 
worked with SATs. This is clearly the recommended 



strategy to follow: implementing the automatic gener- 
ation of the SATs is relatively easy (by systematic ap- 
plication of the RS symmetries as outlined above) and 
results in more transparent and easier-to-construct mod- 
els. Thus, we will refer to SATs when describing the 
effective models for ABO3 compounds; the relevant ones 
(i.e., their representative couplings) are listed in Tables U 
and HH 



4- Long-range interactions in insulators 

The potentials described above can in principle involve 
interatomic interactions of arbitrary spatial range. How- 
ever, in practice we will truncate the spatial extent of 
such interactions, which will constitute one of the approx- 
imations in our models. Generally speaking, such a trun- 
cation can be expected to work well in metals, where the 
free charges provide an efficient means of screening. In 
contrast, the truncation is not justified when we deal with 
semiconductors or insulators, where long-range (strictly 
speaking, infinite-range) Coulomb interactions must nec- 
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TABLE I. Representatives of the symmetry-adapted terms 
(SATs) that couple first-nearest neighbors in the ABO3 cu- 
bic structure. The atom labels correspond to those in Fig. [3] 
(Note that all these representative couplings can be chosen 
so that the two atoms involved are in the same crystal cell.) 
For the atomic displacements we use the compact notation de- 
scribed in the text. We number the couplings to refer to them 
easily in the text. This also allows for a compact notation for 
the coupling coefficients; for example, the SATs sketched in 
Fig.0correspond to coefficients -K#is (left) and if #11 (right). 
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TABLE II. Same as in Table |TJ but involving couplings be- 
tween strains (linear) and atomic displacements of nearest- 
neighboring atoms (quadratic). Strains are given in the Carte- 
sian notation rj a to facilitate the interpretation of the terms. 
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essarily be considered 41 Fortunately, such couplings have 
a well-known analytic form in the limit of long distances, 
and they can be conveniently treated in a way that is 
essentially exact. 

To understand the role of ion-ion Coulomb interactions 
in insulators, let us consider two separate effects. (In the 
following we will implicitly consider the case of short- 
circuit boundary conditions, which corresponds to the 
ideal situation for an infinite bulk material. The treat- 
ment of different electrostatic boundary conditions is dis- 
cussed e.g. inRef.HU.) First, these interactions give raise 
to the so-called Madelung field that contributes to deter- 
mine the cohesive energy of the material. In our model 
potentials, such a Madelung field is captured in the en- 
ergy of the RS. Hence, by taking i?Rs directly as a result 
of the first-principles calculations, we avoid the need to 
model the Madelung energy, as well as the other effects 
(e.g., atomic and short-range interactions associated with 
chemical bonding) that control the basic cohesive energy. 
Second, the Coulombic interaction between ions also in- 
fluences the energy changes associated with the distor- 
tions of the RS. To leading order in the Taylor series, 
such an effect is essentially captured by the electrostatic 



interaction between the dipoles that appear when ions 
move from their RS positions. Such atomic dipoles are 
usually written, within a linear approximation, as 

d ia = Z hcc U i& i ( 21 ) 

P 

where Z* is the so-called Born effective- charge tensor or 
dynamical- charge tensor for atom i. (Strictly speaking, 
we should talk about dipole differences. Yet, here we 
will assume that these local dipoles are zero in the RS, 
which will be the natural choice in most cases.) Note 
that the Born charge Z% a quantifies the dipole caused 
by the displacement of the ionic charge associated with 
ion i at its RS position, as well as all the additional effect 
arising from the electronic rearrangement that occurs in 
response to the atomic distortion. In the case of insulat- 
ing ABO3 perovskites like PbTiC>3 and SrTiC>3, the elec- 
tronic effects are very large and result in Born charges 
that even double the value corresponding to the rigid-ion 
limits 3 - Such huge dynamical charges reflect changes in 
the oxygen-cation bonding that play a crucial role in the 
ferroelectric and response properties of those materials^ 
Hence, when working with insulators, it will be con- 
venient to split the energy terms involving atomic dis- 
tortions Ui into short-range ("sr") and long-range ("lr") 
parts. Thus, for example, we have 

K {n) = K (n),sr + R (n)M ^ 
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for the couplings in E p , where it is important to note 
that the decomposition can be done at all orders in the 
Taylor series. Analogously, the strain-phonon terms in 
E sp can be split as 

Of course, analogous splittings can be considered for the 
parameters that appear in our displacement-difference 
representation of Eqs. (fTU]) and 

Here we will only discuss the lowest-order dipole-dipolc 
interactions, which are captured by the harmonic cou- 
plings K^ 1 and have been described in detail in the lit- 
erature. (Harmonic couplings involving other terms in 
the multipole expansion of the electrostatic energy also 
exist; as usually done in first-principles treatments, we 
will neglect their contribution to the long-range part of 
the energy, and effectively capture their possible effects 
in the short-range part.) Following Gonze and Lee,— we 
write the long-range couplings as 



where 

A Q = ^(e-^Arv, (25) 

& 

and 

D = V A ■ Ar , (26) 

with 

Ar = Rj + tj - Ri - t, . (27) 

This is the usual expression for the Coulombic interac- 
tion between two dipoles, generalized for a medium that 
presents an arbitrary dielectric tensor e M quantifying the 
purely electronic (frozen-ion) response of the material. 
As discussed by Gonze and Lee^ Eq. (j!Mf captures the 
non-analytical behavior of the phonon bands for homoge- 
neous (q = 0) distortions, and the related electrostatic ef- 
fects (e.g., the so-called longitudinal-optical-transversal- 
optical splitting of the phonon frequencies). It is also 
trivial to show that the ASR for the K^' lr coefficients 
translates into the condition 

Z* a p = NJ2 K a p =0 Va, , (28) 

i k 

which guarantees that no net dipole is created by a rigid 
displacement of all the atoms in the crystal. Note that 
the Z* tensors are cell- independent, which allows us to 
use the notation Z* K . Finally, let us mention that in 
an actual atomistic simulation, which usually involves 
a periodically-repeated simulation box or supercell, such 



infinitely-ranged couplings can be accurately computed 
by performing an Ewald summation, as described e.g. in 
Ref. OS 

In this work we only considered the Coulombic dipole- 
dipole term associated with K^K Indeed, as discussed 
below, higher-order long-ranged couplings in E p , and fur- 
ther interactions involving strain in E sp , were either ne- 
glected or treated in an effective way. These approxi- 
mations, which follow the spirit of the usual effective- 
Hamiltonian approach to perovskitc oxides, will be dis- 
cussed in Sections [II C 31 and IIII1 

5. Miscellaneous remarks 

We conclude this Section by commenting on various 
aspects of the model potentials just described. 

Approximations involved - Typically to construct an 
effective potential for a material, one starts by consid- 
ering the simplest possible model that makes physical 
sense, and then extends it only as much as needed to get 
a sufficiently accurate description of the first-principles 
data of interest (i.e., a good description of what is usually 
called the training set of first-principles results). Given 
the conceptual simplicity of our proposed potentials, it 
is straightforward to identify three qualitatively differ- 
ent ways in which they can be systematically extended. 
Indeed, our models can be improved as regards (i) the or- 
der of the polynomial expansion, (if) the spatial range of 
the interatomic couplings considered, and (Hi) the com- 
plexity of the coupling terms, i.e., the maximum number 
of atoms (bodies) involved in the couplings. These three 
truncations constitute the approximations of our models. 

Relation with effective-Hamiltonian work- For the 
most part, the connections between our method and 
the above-mentioned effective-Hamiltonian approach are 
rather obvious. Yet, there are a couple of subtle points 
that deserve a comment. 

The effective Hamiltonians often include local variables 
that account for the inhomogeneous strains that may oc- 
cur in the material; further, the energy landscape for such 
local strains is typically derived from the elastic constants 
associated with the homogeneous ones, following the ap- 
proximation proposed by Keating^ In our models, inho- 
mogeneous strains are naturally captured by the appro- 
priate atomic distortions {ui}, as illustrated in Figs.[2][c) 
and[^d). The energy changes associated with such local 
strains are given by E pi and there is no need to derive 
them from the elastic constants for homogeneous cell de- 
formations. (Of course, one should note that the force 
constants and the clastic constants are con- 

nected by well-known relations, and the latter can be 
computed from knowledge of the former^ 7 - As explained 
in Section III CI we include in our models the exact first- 
principles results for both and C^ 2 \ so that the 
relations between such coefficients arc fulfilled by con- 
struction.) Additionally, our models also capture cor- 
rectly the energy changes associated with strong strain 
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gradients. This is a definite improvement over the usual 
effective-Hamiltonian approach, especially when taking 
into account the growing interest in flexoelectric effects 
resulting from large strain gradients near ferroelectric do- 
main walls j2I etc. 

Secondly, the action of an external electric field £ can 
be trivially incorporated in an effective-Hamiltonian sim- 
ulation by including the leading coupling term between 
the field and the local dipoles that are the basic variables 
of the models Equivalently, within our approach (and 
as long as we are dealing with insulators) , we can use the 
effective-charge tensors Z* to compute the local dipole 
di resulting from atomic displacements [Eq. ([2~Tj)], and 
write the corresponding energy as 



E eS ({ui},ri;£) = E cS ({ui},r)) - 2J d ia £ a 



(29) 



Finally, let us note that the action of an external stress 
or pressure <x can be treated in an analogous way^ by 
introducing 

£ e ff(W,»7;<7) = E eS ({u i },r l )+N^2r) a a a . (30) 

a 

Here the sign convention is chosen so that a positive stress 
implies a compression of the material. 

Implementation in a simulation code .- Let us briefly 
mention some details of our implementation of a sta- 
tistical simulation [in particular, a Monte Carlo (MC) 
scheme] based on our model potentials. 

First, let us note that the long-range part of the force 
constant matrix, K^' lr , depends on the specific size and 
shape of the periodically-repeated simulation box used 
for the MC runs. Hence, for a given supercell, we com- 
pute these parameters before the MC simulation starts, 
by performing the corresponding Ewald sums that take 
into account interactions between periodically-repeated 
dipole images. Then, we add up the long-range and 
short-range parts of K^ 2%> to obtain a total harmonic in- 
teraction term that effectively couples all atoms in the 
simulation box. This is what we use for the energy eval- 
uations in the simulation. 

Once we have a supercell-depcndent potential, the un- 
derlying lattice-translational symmetry allows us to store 
only the interactions between the atoms in one elemental 
unit cell and all other atoms in the supercell. Hence, the 
storage requirements grow linearly with number of unit 
cells in the supercell. 

In our MC simulations, we attempt to change the 
strains only after completing one sweep through all the 
atoms in the simulation supercell. It is therefore conve- 
nient to recalculate the parameters controlling the ener- 
getics of the atomic displacements, such as for example 



K 



(2) 

iotj/3 



(31) 



after the strains are updated. These strain-dependent 
parameters are then used for energy evaluations during 
the sweep over atomic displacements. 



The SATs for the calculation of the anharmonic part 
of our models are automatically generated based on the 
symmetry of the RS. We store them in symbolic form, so 
that they can be used both for the calculation of energy 
and (via a simple manipulation of the polynomial) forces 
on the atoms. 



C. Parameter Calculation 

Once we have defined a potential, many schemes can 
be applied to calculate its parameters. Here we describe 
the strategy that we followed in this first application of 
our effective models, which takes advantage of the di- 
rect availability of first-principles results for many of the 
terms in the potential. Some approximations that we 
used for the treatment of long-range interactions, which 
are somewhat specific to the case of insulators undergo- 
ing structural phase transitions, are also described. 



1. Parameters computed directly from first principles 

The low-order couplings of our model potentials quan- 
tify the response of the RS of the material to small pertur- 
bations, may they be atomic distortions, cell strains, or a 
combination of both. In particular, the leading harmonic 
terms K^ 2 \ C^ 2 \ and A^ 1,1 ' can be obtained directly 
from density-functional perturbation theory (DFPT) cal- 
culations as those described in Refs.[35, 39 and 30|. DFPT 
schemes are efficiently implemented in widely available 
first-principles codes, such as the Abinit package^ used 
in this work. Alternatively, one could obtain the same 
information by performing systematic finite-difference 
calculations considering both atomic displacements and 
strains. Such an approach, which is somewhat more ele- 
mentary but equally valid, is available in all major first- 
principles packages. Hence, we can conclude that com- 
puting exactly the harmonic parameters of models like 
ours is a trivial task nowadays. 

Let us stress that the ability to incorporate an exact 
description of the harmonic energy of the material by con- 
struction is a great asset of our models. Indeed, in most 
materials the thermodynamic properties are essentially 
captured at the harmonic level, with small corrections 
coming from anharmonic effects; hence, a good descrip- 
tion of the harmonic lattice-dynamical properties is crit- 
ical. Further, even in cases with soft-mode-driven phase 
transitions, it is the harmonic part of the energy what es- 
sentially determines the nature of the leading structural 
instabilities. Hence, also in such situations, a faithful 
harmonic description seems mandatory to have an ac- 
curate model. Figure 2] shows representative results for 
our model of PbTi03. As we can see, the description 
of the force-constant bands of the RS is exact, and the 
small discrepancies between the shown density-of-states 
(DOS) plots come from differences in the way BZ inte- 
grations are performed in Abinit and in our codes. The 
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DOS (LDA) DOS (model) 




FIG. 4. (Color on line.) Left: Dispersion bands of cubic 
PbTiOs, as calculated from first principles (lines) and ob- 
tained from our effective model (circles). The bands corre- 
spond to the eigenvalues X q j of the Fourier-transformed force- 



constant matrices K 



(2) 



which we call stiffness coefficients. 



The leading structural instabilities are labeled (see the text), 
and sketched in Fig. [5] The color code indicates the dominant 

( 2) 

atomic character of the K q eigenvectors. Right: Density of 
states (DOS) plots constructed from the Kg eigenvalues, as 
obtained from first-principles simulations using a very fine q 
point mesh, and from our effective potential by solving the 
eigenmode problem for an 8x8x8 supercell and making use 
of a simple interpolation between the computed eigenvalues. 



most important structural instabilities, marked in Fig. |4] 
and sketched in Fig. [5l are also reproduced exactly. 

As regards the anharmonic terms, one could try a sim- 
ilar direct calculation of each one of the parameters. For 
example, to compute the strain-phonon couplings A*- 1 ' 2 -*, 
one could run DFPT calculations for the RS subject to 
a small strain 5r). The resulting force-constant matrix 
would be described in our model by 



K 



(2) 
iajP 



which would allow us to calculate the targeted couplings. 
Following a similar scheme - e.g., by running DFPT 
calculations of distorted configurations in which some 
atomic displacements are frozen in - one could access 
the parameters in -B a nh- 

As described below, we tried such an approach when 
constructing our models for PbTiC>3 and SrTi03, specifi- 
cally in what regards the strain-phonon couplings. Based 
on our experience, we believe that such a systematic 
scheme may render accurate potentials in relatively sim- 
ple cases, i.e., whenever the RS does not present struc- 
tural instabilities. On the other hand, in the challenging 
situations here considered, this strategy may be imprac- 
tical if a very precise description of some PES features 
is targeted. Indeed, we found that the PES of materials 
like PbTiO,3 or SrTiC>3 is strongly anharmonic; more pre- 
cisely, if we aimed at an accurate description of the whole 
PES connecting the RS with the lower-energy phases, we 
would need to consider a Taylor series extending up to 
a rather high order. In such cases it seems more con- 
venient to adopt an effective approach, aiming at rcpro- 





(a) FE 2 



(b)AFD^ 



(c)AFD^ 



FIG. 5. Sketch of the atomic displacements corresponding 
to the most important structural instabilities in ABO3 per- 
ovskite oxides. Panel (a): Ferroelectric instability. Panel (b): 
Anti-ferrodistortive instability with neighboring Oe octahedra 
along the z direction rotating in phase. Panel (c): Same as 
in (b), but with octahedral rotations modulated in anti-phase 
along z. 



ducing the PES only around the RS and the most rele- 
vant low-energy structures. This permits a lower-order 
expansion that quantitatively captures the main effects 
and retains much of the physical transparency of the sim- 
pler (effective-Hamiltonian and phenomcnological) mod- 
els traditionally used to investigate phase transitions, 
which include only as many terms as strictly needed for 
a qualitatively correct description. 



2. Parameters fitted to first-principles results 

To compute the higher-order couplings of our effee- 

tive potentials - i.e., K with n > 2 and A with 
m + n > 2 -, it is convenient to implement a fitting 
procedure aimed at obtaining a model that reproduces a 
training set of first-principles results. Here we describe 
the strategy we adopted in our work with PbTi03 and 
SrTi03, where the training set was composed of low- 
energy structures that are more stable than the RS, and 
the key properties that we request our models to cap- 
ture are energy differences and equilibrium atomic con- 
figurations. Nevertheless, the ideas presented are rather 
general and can be easily adapted to other situations. 

In essence, our parameter-optimization calculations 
were based on three goal functions defined in the fol- 
lowing way. Let the superindex s number the structures 
({uf }, T7 S ) in our training set. First, to get our model to 
reproduce the first-principles energies {E s }, we consid- 
ered the goal function 



[E cft [P]({ut},V s )-E s 



(33) 



where V represents all the free adjustable coefficients in 
the model and the parametric dependence of -E e ff on V is 
indicated. Second, all the structures in our training sets 
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were stationary points of the PES (minima or saddles). 
Hence, we imposed the zero-gradient condition for such 
structures by minimizing the goal function 

GFveP) = E l|V£ off ™<},rf )|| 2 , (34) 

s 

where the gradient includes derivatives with respect to 
both atomic distortions and cell strains. Finally, aim- 
ing at an improved description of the lattice-dynamical 
properties of key low-energy structures, we also used a 
goal function that contains information about the corre- 
sponding Hessian matrices. More precisely, we used 

OF hcss (7>) = E E P'Plte). (35) 

where {q} s is a set of g-points of the first Brillouin zone 
of structure s (we restricted ourselves to zone-center and 
zone-boundary g-points). The function T> s [V](q) quan- 
tifies the difference between the Hessian for structure s 
obtained from the model (K^g q ) and its first-principles 
counterpart (K s q ); we define it as 

V s [V]{q) = E \\K s cS . q v s q] - A^.f , (36) 

3 

where vtj and X q j stand, respectively, for the eigenvec- 
tors and eigenvalues of the first-principles Hessian K q . 
This strategy to compare the Hessian matrices allowed 
us to achieve meaningful parameters in a reliable and ro- 
bust way; in contrast, we found that simpler schemes, 
based on a direct comparison of eigenvalues or eigenvec- 
tors, lead to difficult optimization problems that present 
many spurious local minima of the goal function. 

The above functions can be combined to run opti- 
mizations targeting simultaneously at different proper- 
ties. However, it is not clear a priori how to weight 
the different goal functions in order to construct a sin- 
gle QF that renders a well-posed optimization problem. 
Hence, we adopted the following alternative approach, 
which we used to generate most of the results presented 
in Section IIIII We start the parameter optimization by 
minimizing one of the goal functions, QF\. Then a sec- 
ond goal function QF2 is minimized, with the parameters 
subject to the constraint that the result for QF\ must be 
preserved within a certain tolerance. In this way, suc- 
cessive optimizations can be performed, with constraints 
involving all previously-optimized goal functions, until 
we impose all the necessary conditions. Naturally, the 
tolerances for the constraints can be chosen so that the 
most critical properties are reproduced better. Typically, 
in our work with PbTi03 and SrTiC>3 we started by min- 
imizing QF Ei as we prioritize that our models reproduce 
correctly the first-principles energies of the structures in 
the training set. Then, the most usual sequence of opti- 
mizations involved QF^e, <?-Fhess evaluated at the T (i.e., 
q = 0) point of the lowest-energy structure(s), and finally 
QFhoss evaluated at selected zone-boundary g-points of 
the lowest-energy structure(s). 



The optimization of GFh css was never prioritized in the 
applications considered in this work; in fact, we found 
that, when working with relatively simple (low-order) 
models as the ones considered here, it is not realistic 
to aim at a very precise description of the first-principles 
Hessians of structures that deviate significantly from the 
RS. Nevertheless, we found that it was often possible 
to adjust the low-lying eigenmodes at a reduced num- 
ber of g-points. Also, generally speaking, we found that 
considering GFh css was a good strategy to obtain energy- 
bounded potentials, as such an optimization step helps 
to impose the stability of the ground state structure. 



3. Further comments on the long-range interactions 

As mentioned above, the atomic interactions in insula- 
tors can be conveniently decomposed in short- and long- 
range parts. Further, at the harmonic level we have a 
simple analytical expression for the dipole-dipole cou- 
pling [Eq. (|24p] that depends on the RS geometry, the 
dynamical charges Z* , and the dielectric tensor e x . Con- 
veniently, these tensors, as well as the the decomposition 
of into K^' sr and K^' lr , are produced automati- 
cally by most DFPT implementations; in particular, they 
are readily provided by Abinit. [The typical DFPT 
scheme computes the total interatomic force constants. 
Then, in essence, it is assumed that the long-range part 
K^ ,lr is given by the dipole-dipole term in Eq. (|24p . 
and the short-range part is obtained as K^' sr = 
— K^' lr .] Alternatively, all the relevant parameters con- 
trolling the dipole-dipole interactions can be obtained by 
considering the response to finite electric fields4i 

As regards the anharmonic terms, we could continue 
to distinguish between short- and long-range couplings. 
In essence, the anharmonic long-range couplings in E p 
would capture the changes in the effective charges or 
dielectric constants that may be caused by the atomic 
displacements and which affect the magnitude of the 
Coulombic dipole-dipole interactions. As regards the 
strain-phonon couplings in E sp , an additional effect 
comes from the change in the cell shape and dimensions. 

Our model potentials provide a framework to capture 
such effects by considering appropriate high-order terms. 
Unfortunately, considering such couplings would result in 
computationally-heavy atomistic simulations. Indeed, as 
discussed in Section lllB 5[ for a practical implementation 
of the harmonic long-range interactions it is convenient to 
precalculate, for the RS geometry and our specific choice 
of simulation supercell, the dipole-dipole couplings by 
performing the appropriate Ewald sums. Once the in- 
teraction coefficients K^' lr are known, the correspond- 
ing energy can be readily obtained during the course of 
the simulation; yet, because such a term couples all the 
atoms in the supercell, its calculation is by far the most 
time-consuming part of the energy evaluation. In prin- 
ciple, one may proceed similarly with the higher-order 
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long-ranged terms. For example, consider 

A (l,a) =A (l,2),«- + A (l I 2),Jr j (37) 

which is the leading strain-phonon coupling for mate- 
rials like PbTi03 and SrTi03. In this case, we can 
model A*- 1 ' 2 -''''" by considering the dependence of K^ ,lr 
[Eq. ([24)) ] on a strain r\ to linear order. (To do this, one 
could proceed by introducing the strain dependence of 
the effective charges, dielectric tensor, and interatomic 
distances in Eq. (f2~4")l . and then Taylor expand with re- 
spect to rj.) The corresponding coefficients could be pre- 
computed for the RS geometry and particular simulation 
supercell, which would permit an easy (but still computa- 
tionally costly) evaluation of such an energy contribution 
during the course of the simulation. 

In our work with PbTiC>3 and SrTiC>3, wanting to ob- 
tain models that allow for fast simulations, we did not 
treat explicitly the anharmonic corrections to the long- 
range dipole-dipole interactions. Yet, we captured the 
effects on the properties of interest (e.g., the energy, 
equilibrium structure, and Hessian of low-lying phases) 
in the short-range anharmonic couplings. Whenever the 
anharmonic couplings are determined by the fitting pro- 
cedure outlined in Section III C 21 this can be done in the 

— -(n) ,lr 

most natural way. We simply assume that K = 

~ (n,m),lr _ - — (n),sr 

A =0, and fit the anharmonic terms K and 

— (n,m),sr 

A to reproduce first-principles information about 

the structures in our training set, thus capturing effec- 
tively the consequences of possible anharmonicities in the 
long-range couplings. 

Additionally, we also computed the strain-phonon cou- 
plings directly, without performing any fit, by proceeding 
in the following way. We considered the full interatomic 
constants for the RS and strained configurations, and as- 
sumed that the following approximate version of Eq. (|32[) 



K 



Si] 



(2),Zr 
iaj/3 



5rj 



VA (1 ' 2) ' sr fe 



(38) 



holds within a certain spatial range (i.e., for a maxi- 
mum separation of atoms i and j). Then, we demanded 
that the short-range part of A^ 1 ' 2 ' capture strain-induced 
changes in both K (2) ' sr and K^' lr . It must be noted 
that, because of the spatial truncation, the A.^ 1 ^' sr 
thus calculated will in general break translational invari- 

ance. To remedy this, we added to A^ 1 ' 2 -*' 51 * a correction 
AA (i,2), S r that 

was determined by demanding that our 
model reproduce exactly the Hessian of the strained con- 
figurations at the r point. In this way, by imposing a 
correct description of the acoustic modes, we restore the 
ASR. Further, this procedure also guarantees that the 
effect of strain on the T distortions, which are critical 
for the investigation of ferroic perovskites like ours, is 
captured by our models. As shown in Section IIIIB1 this 



approximation leads to a very precise description of the 
strain effects on the force-constant bands in the case of 
PbTi0 3 . 



III. EXAMPLES OF APPLICATION 

Now we describe the model potentials for ferroic per- 
ovskites PbTi0 3 (PTO) and SrTi0 3 (STO) that we con- 
structed following the above scheme. These materials are 
representative of the large family of compounds under- 
going structural phase transitions driven by soft phonon 
modes. The lattice dynamical properties of such sys- 
tems are strongly anharmonic, and the description of 
their transitions requires the use of high-order potentials. 
Further, in the case of these perovskite oxides the rele- 
vant energy scale for the soft mode instabilities is rela- 
tively small, of about 50 meV per formula unit (f.u.) or 
less. Hence, achieving a good description of such com- 
pounds constitutes a challenge for first-principles theory 
and, naturally, for our model-potential approach. 

Additionally, PTO and STO present peculiarities that 
make them especially interesting in the present context. 
At temperature Tq = 760 K, PTO undergoes a transition 
between the high-T paraelectric structure (i.e., the ideal 
cubic perovskite prototype, with space group Pm3m, 
that we take as our RS) and its low-T ferroelectric (FE) 
phase (with tetragonal space group P4wm) i 21 ' 42 The 
structural distortion that appears at low temperatures 
has a polar character, and it essentially involves a dis- 
placement of the Ti and Pb cation sublattices against the 
06-octahedron network, as sketched in Fig. Eta). Note 
that this corresponds to the condensation of a soft mode 
at the zone center (at the T point) of the BZ of the RS. 
This transition has a significant first-order character that 
previous theoretical work has linked with the accompa- 
nying deformation of the cellfii further, first-principles 
theory predicts that cell strains are critical to determine 
the symmetry of the ground state of PT042 Hence, to 
model this compound we have to deal with both the FE 
instability responsible for the transformation at Tq and 
the strain-phonon couplings that have a strong impact in 
the occurring equilibrium phases and the features of the 
FE transition. 

SrTi03 too undergoes a single phase transition, as it 
transforms at 105 K from the high-T cubic perovskite 
phase to a low-T structure of tetragonal (74 /mem) 
symmetry^ The structural distortion occurring in the 
low-T phase involves concerted rotations of the Og oc- 
tahedra about the tetragonal axis, with the peculiarity 
that Og groups that are first neighbors along z rotate in 
antiphase. Such a pattern is denoted a°a°c~ in the well- 
known notation introduced by Glazer^i and corresponds 
to a so-called antiferrodistortive (AFD) mode associated 
with the R point of the BZ of the RS [q R = 7r/eto(l, 1, 1), 
where do is the lattice constant of the RS cubic unit cell] ; 
the corresponding atomic displacements are sketched in 
Fig. EJc). Additionally, STO is close to presenting a 
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FE instability; in fact, this compound is experimentally 
believed to be a quantum-paraelectric, i.e., a material 
whose ferroelectricity is suppressed by quantum fluctu- 
ations (i.e., the wave-like character) of its constituting 
atoms. 14 i 45 i 46 Further, previous first-principles work has 
shown that the FE and AFD soft modes compete in 
STOji 3 - complicating even more the description of the 
behavior of the material at low temperatures. Hence, 
STO offered us the possibility of testing our approach in 
cases in which several structural instabilities are relevant 
and their interaction must be considered in detail. 

We first describe our work with PTO, which turned out 
to present all the challenging features that we had antici- 
pated (i.e., the very critical strain-phonon couplings) and 
additional ones that we were not expecting (i.e., a very 
significant competition between FE and AFD modes). 
Hence, we discuss the case of PTO in detail, giving il- 
lustrative examples of how our models can be extended 
when it is necessary to do so. In contrast, it was rel- 
atively easy to obtain a sound model for STO. Hence, 
in that case we will present a very minimal approach to 
the construction of an effective potential. In both cases, 
we will describe the T-driven transitions obtained when 
solving our models by means of Monte Carlo simulations 
(in which, as usually done, we treated atoms as clas- 
sical objects), showing that they capture correctly the 
basic experimental behaviors. We will also comment on 
the probable origin of the quantitative discrepancies ob- 
served between our model predictions and experiment. 
Note that here we will not elaborate much on the physics 
emerging from our models, as such a discussion falls be- 
yond of the scope of this paper. 



A. First-principles and Monte Carlo methods 

All first-principles calculations were done with the 
Abinit package^ and employed the local-density ap- 
proximation (LDA) to density functional theory 47 ' 48 
The ionic cores were treated by using extended norm- 
conserving Teter pscudopotcntials^ and the following 
electrons were considered explicitly in the calculations: 
Pb's 5d 10 , 6s 2 , and 6p 2 ; Sr's 4s 2 , 4p 6 , and 5s 2 ; Ti's 3s 2 , 
3p 6 , 3d 2 , and 4s 2 ; and O's 2s 2 and 2p A . Electronic wave 
functions were represented in a plane-wave basis trun- 
cated at 1500 eV. We used an 8 x 8 x 8 fc-point grid 
to compute integrals in the Brillouin-zone of the 5-atom 
perovskite cell, and equivalent meshes for other cells. In 
structural relaxations, atomic positions were optimized 
until residual forces on atoms were below 10 -4 eV/A. 
The interatomic force constants, elastic constants, Born 
charges, and dielectric tensor were calculated by using 
the DFPT implementation in Abinit. The K s q matri- 
ces, from which the real-space interatomic constants are 
obtained, were computed for a 2 x 2 x 2 q-point mesh; 
in agreement with previous studies^ this was found to 
be sufficient to get accurate results. The resulting cut-off 
radius for the short-range interactions is therefore about 



6.8 A for all of the presented models. 

Thermal averages of the quantities of interest were cal- 
culated by a standard Metropolis Monte Carlo method* 5 - 1 - 
The Markov chain was constructed by sequentially con- 
sidering movements of (i.e., by sweeping through) all 
atoms in the simulation box. After each sweep, a sin- 
gle attempt to modify each of the strain components 
was made. Both the attempted displacements and strain 
modifications were drawn from appropriate uniform dis- 
tributions, whose widths were varied by a simple linear 
controller with the goal of attaining an acceptance rate 
of 50% on average. In most cases, we used an 8 x 8 x 8 
periodically-repeated simulation supercell, and thermal- 
ized the material by running 20000 MC sweeps starting 
from the RS (i.e., Ui = rf = 0). The averages for the 
relevant structural distortions were then calculated from 
20000-40000 additional sweeps, and we checked conver- 
gence by inspection of the corresponding histograms. At 
temperatures in the vicinity of the phase transitions, this 
procedure did not lead to converged results because of 
either slow thermalization or finite-size effects. In such 
cases, we found it necessary to run the calculations of up 
to 80000 sweeps in 10 x 10 x 10 simulation supercells. 
For presentation purposes, all of the computed average 
distortions were rotated so that the axes for the FE po- 
larization and AFD rotations lie along the [001] (resp. 
[Ill]) Cartesian direction for tetragonal (resp. rhombo- 
hedral) phases. 

B. PbTi0 3 

1. Harmonic terms i?har({iii}) and E s (ri) 

The first step in the construction of our model poten- 
tial is the computation of the harmonic energy terms, 
i?har and E s , for which we use the DFPT schem o 30 ' 35 im- 
plemented in Abinit4S (As mentioned above, the third 
harmonic term - i.e., the strain-phonon coupling A^ 1,1 -* 
in E sp - is identically zero in PTO and STO due to the 
cubic symmetry of the RS.) Representative results are 
given in Fig. 01 which shows the bands corresponding to 
the stiffness coefficients or force constants of the cubic 
RS. (These are the eigenvalues X q j of the Hessian matri- 
ces K c ff.q introduced in Section lll C 21 ) Notably, we find 
that some distortions have a negative stiffness, indicat- 
ing that they are structural instabilities of the RS. The 
leading instabilities arc pictorially represented in Fig. [5] 
the FE soft mode at the T point [panel (a)], the in-phase 
AFD mode at the M point [q M = 7r/ao(l, 1,0), where ao 
is the lattice constant of the RS cubic unit cell; panel (b)], 
and the antiphase AFD mode at the R point [panel (c)] . 
As shown in Fig. 21 the energetics of all such instabilities 
is captured exactly, at the harmonic level, by our model. 

As a result of our DFPT calculations, we obtained 
an £"har term that includes all short-range interactions 
within a spatial range slightly below 7 A. (For ex- 
ample, this includes couplings between Ti pairs that 



16 



12 
10 



•< 
> 

& 6 

< 



~i 1 1 r 



FE Z + x □ 
FE >yz + x □ 
AFD' Z + x □ 
Pb TiO 



3 4 5 6 7 8 
Interatomic distance (A) 



FIG. 6. (Color on line.) Deviation of the interatomic force 
constants calculated for selected PTO distorted structures 
from the RS results. Each point quantifies the difference be- 
tween the 3x3 force-constant matrices, for a specific atom 
pair, computed for the reference and distorted structures. The 
results are shown as a function of interatomic distance; for 
clarity reasons, the interatomic distances are shifted slightly 
to reduce overlap. Note that, when two different atoms are 
involved in the pair, we overlap the corresponding symbols; 
thus, for example, crossed squares correspond to pairs involv- 
ing Ti (cross) and O (square) atoms. 



are 3rd nearest neighbors.) Additionally, £har includes 
the already-mentioned analytic form of the long-range 
dipole-dipole couplings, which involves 5 symmetry- 
independent parameters [i.e., 4 Born effective charges 
(which reduce to 3 independent ones if the ASR in 
Eq. (|28|) is considered) and 1 dielectric constant that fully 
defines the diagonal and isotropic tensor] . As regards the 
harmonic elastic constants in E s , the model incorporates 
the 3 symmetry-independent terms that define the full 
elastic tensor for a crystal with cubic m3m point sym- 
metry. 



2. Fitting E &ah ({ui}) 

Next we tackled the construction of the anharmonic 
terms of the potential. We first considered the case in 
which the cell is fixed to be that of the RS, i.e., we as- 
sumed r\ = and focused on i? a „h- As described above, 
we computed E an h by fitting its parameters to a set of 
relevant first-principles data. Naturally, we populated 
our training set with information about the low-energy 
structures that can be accessed by condensing the differ- 
ent instabilities of the RS. More specifically, our list of 
low-symmetry phases contains FE structures of tetrag- 
onal (FE Z ) and rhombohedral (FE xyz ) symmetries, as 
well as several AFD-distorted phases (AFD°, AFD® yz , 
and AFD*). [We use the notation of Fig. [5l with the 
xyz subscript denoting the simultaneous occurrence of 
a distortion type along/about the three Cartesian axes, 
and with the same amplitude for the three of them.] In 
addition, we also considered a hybrid structure, FE^^ + 




FIG. 7. (Color on line.) Potential-energy wells connecting the 
RS of PTO with the low-symmetry phases defined in the text. 
The results obtained from our model potential are shown with 
lines, and the points indicate the first-principles results for the 
energy minima or saddles. All the states shown preserve the 
cubic cell of the RS (r) — 0). The amplitudes |u| in Angstrom 
correspond to collective distortions involving several atoms. 
The AFD" and AFD^;, curves are essentially on top of each 
other and cannot be distinguished. 



AFD" y2 , that has rhombohedral i?3c symmetry and in 
which both FE and AFD patterns occur simultaneously. 
Note that we determined such low-symmetry structures 
ab initio by (i) distorting the RS according to a spe- 
cific soft-mode eigenvector and then (ii) using this as the 
starting point of a structural relaxation that preserves 
the symmetry of the initial configuration. In all cases we 
computed the equilibrium structure, energy, and Hessian 
matrix. The relevant structural parameters of the con- 
sidered phases are given in Table HEH together with the 
energies relative to the RS. 

In order to fit -E an hj we worked with the displacement- 

~(n) 

difference representation and K parameters (with n > 
2) of Eq. (fTU|) . We restricted ourselves to models that in- 
clude only pairwise interactions and extend up to 4th 
order in the Taylor series. These approximations define 
the minimal model needed to capture structural phase 
transitions like the ones we want to describe, and are 
analogous to the ones adopted in most of the previous 
theoretical works that we are aware of. (One of the few 
exceptions is the inclusion of high-order terms for the 
local polar modes considered in Ref. [ill) In our case, 
we maintained such approximations in order to keep our 
models relatively simple and computationally efficient, 
as well as to test the actual ability of such an elemen- 
tary potential to reproduce the first-principles data in a 
quantitative way. 

As regards the spatial extent of the anharmonic cou- 
plings, most of the previous works on phenomenological 
models and effective Hamiltonians adopt what is some- 
times called the on-site anharmonicity approximation, 
which implies that the non-harmonic couplings are taken 
to be strictly confined in space and contribute only to 
the self-energy of the atoms or local modes iS r 11 ' 20 ' 21 In- 
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TABLE III. Structural parameters of the considered low-symmetry PbTiOs structures (see the text) calculated for a cubic 
cell with lattice constant a = 3.880A. First-principles LDA results are presented along with those obtained from the model 
discussed in the text. We show the T-point displacements corresponding to the polar distortion, vf Ka , given in Angstrom; we 
chose them so that there is no rigid shift of the whole structure (i.e., upbz + «Tiz + 2uq\ z + uo3z = 0). For P4mm, all the 
displacements are along the z direction and we have Uq 1z = Uq 2z - For Rim and R3c, we have u^ x — u^ y — u^ z for the Pb and 

Ti atoms, as well as Wo3z = «oii = u 02 y an d Uoiz — u oi y — u 02x = u 02z = u 03x = u 03y The amplitude of the AFD modes 
is quantified by the corresponding Oe rotation angle given in degrees (in the rhombohedral cases, we have equal-magnitude 
rotations about the three Cartesian axes; we give the rotation angle for one axis.). Energies are given in meV/f.u., taking the 
result for the RS as the zero of energy. 
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tcrcstingly, our first-principles results give us a direct way 
to test whether such an approximation is justified. Fig- 
ure IH1 shows the difference between the harmonic inter- 
atomic couplings computed for the RS (which are given 
by K 1 - 2 ^ directly) and those corresponding to the sev- 
eral distorted states of PTO that maintain the cubic cell 
(which are described by plus a distortion-dependent 

~(n) 

correction involving K with n > 2). From these re- 
sults, it is apparent that the distortion-induced changes 
decay very rapidly with the interatomic distance, indicat- 
ing that the anharmonic corrections have a limited spa- 
tial range; similar calculations for other distorted config- 
urations confirmed this conclusion. Hence, our model for 

~(«) 

PTO included only anharmonic K couplings between 
neighboring atom pairs (i.e., each Pb atom is coupled 
with its 12 neighboring oxygens, and each Ti atom with 
the 6 oxygens in the surrounding Og group), which re- 
sults in couplings extending up to about 3 A. Note that 
this approximation is essentially equivalent to the on-site- 
anharmonicity assumption of the effective-Hamiltonian 
method, but adapted to our displacement-difference rep- 
resentation. Together with the other truncations men- 
tioned above (pairwise interactions, 4th-ordcr Taylor se- 
ries), this local-anharmonicity approximation results in 
the 15 SATs listed in Table Q] 

Using the model and training set described above, we 
fitted the 15 anharmonic parameters of Table U by suc- 
cessive optimization of the GTe, Q?VEi and QT^ CS& 
goal functions, following the procedure outlined in Sec- 
tion lll C"2l In GFhess we considered the Hessian matrices 
of distorted configurations, including modes correspond- 
ing to the r and, in some cases, R points of the BZ 
of the RS. For each g-point, wc considered only the 6 
lowest-lying optical eigenmodes (i.e., we did not fit to 



the full spectrum). As evidenced by Table Hill and Fig. 
the model thus constructed describes with good accuracy 
our first-principles results for the equilibrium structures 
and energies of the relevant r\ = configurations. Addi- 
tionally, Fig. [5] shows the results that our model gives for 
the force-constant bands of two distorted structures; as 
expected, the low-lying Hessian eigenmodes are reason- 
ably well reproduced, and the inaccuracies grow as we 
move up in energy. 

To test our model for a fixed-cell version of PTO, we 
ran MC simulations and computed the evolution of the 
equilibrium structure as a function of temperature. Fig- 
ure [9] shows our basic results, which reveal a sequence 
of two phase transitions: At T w 200 K the material 
develops a spontaneous polarization, which manifests it- 
self in a non-zero value of the dipole moments averaged 
over all cells in the simulation box. Such a transition 
drives the system from its high-T cubic (Pm3m) phase 
to a rhombohedral (i?3m) one; the spontaneous polar- 
ization is parallel to the rhombohedral axis, which lies 
along the [1,1,1] Cartesian direction. Such a Rim struc- 
ture is usually thought to be the ground state of PTO 
subject to the r\ = condition! 11 i 43 However, our MC 
simulations rendered a second transition, at T w 100 K, 
in which an AFD mode freezes in. More precisely, at low 
temperatures we observe the occurrence of a distortion 
involving antiphase rotations of the Og groups about all 
three Cartesian axes, which we denoted by AFD" yz in the 
description above. The spontaneous polarization remains 
essentially unaltered upon the condensation of this AFD 
mode, and the new phase presents the i?3c rhombohedral 
space group. Remarkably, this structure was not part of 
the initial training set that we used to fit the parameters 
in -Eanh; indeed, we discovered it by running MC simu- 
lations with our initial model potentials for PTO, which 
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(a) cubic (Pm3m) (b) FE Z (P4ram) (c) FE xyz (R3m) 




FIG. 8. (Color on line.) Force-constant bands corresponding to three different PTO structures, all maintaining the cubic RS 
cell. Black solid lines show the results of our model potential, and red dashed lines the first-principles results. 



led us to better characterize it from first principles and 
eventually include it in the training set. This clearly il- 
lustrates the usefulness of our model-potential approach 
to discover new phenomena. [The force-constant bands 
of Fig. |8tc) already indicate that the Rim structure can- 
not be the ground state of our fixed-cell version of PTO. 
Note that the bands for the R3m structure present a neg- 
ative stiffness for some i?-point modes, which correspond 
exactly to the low-T AFD instability observed in the MC 
runs.] 



3. The strain-phonon term E ap ({ui},r)) 

We began by considering in E sp ({ui}, rj) some of the 
lowest-order terms that are not zero by symmetry, i.e., 

f 1 2^ 

those corresponding to the coefficients A k ' ' or, cquiv- 

~(1>2) 

alcntly, A . This constitutes the minimal approxi- 
mation that captures the strain-phonon couplings lead- 
ing to physically relevant phenomena in ferroelectric per- 
ovskites (e.g., piezoelectricity and the elastic effects as- 
sociated with the structural transitions), and is anal- 
ogous to the one adopted in the effective-Hamiltonian 
literature^— 
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FIG. 9. (Color on line.) Temperature-dependent polarization 
and AFD a order parameters of PTO as obtained from MC 
simulations of our model under the rj — condition. 



As regards the spatial extent of the A^ 1,2 -* interatomic 
couplings, the effective-Hamiltonian works have tradi- 
tionally adopted an on-site approximation that is anal- 
ogous to the one used for the anharmonic terms in 
iJ anh jihii consequently, only one-body interactions are 
typically considered. Further, strain effects on the long- 
range dipolc-dipolc interactions have never been treated 
in the literature, to the best of our knowledge. 

In our case, we went beyond such approximations by 
computing A^ 1,2 ^ via the approach described by Eqs. (|32|) 
and (|55|) . using strains of ±2% for the finite-difference 
calculations. (The spatial extent of the A 1 - 1 ' 2 - 1 interac- 
tions thus computed is essentially identical to that of the 
terms.) The model constructed in this way, which 
we call L , captures very accurately the strain depen- 
dence of the force-constant bands, as can be appreciated 
in Fig. HUT a) ; note that such a good agreement through- 
out the BZ validates our approximate method to treat 
the effect of strain on the long-range interactions between 
dipoles, which was described in Section III C 31 Moreover, 
this L° model also renders the correct low-T structure for 
the real (unconstrained) PTO: Indeed, the experimental 
ground state of bulk PTO at ambient pressure is tetrag- 
onal (PAmm space group), as opposed to the rhombohc- 
dral (i?3c) solution that we predict when imposing the 
r\ = condition. Remarkably, the strain-phonon cou- 
plings calculated with our finite-difference scheme cap- 
ture such an effect, even though they were not explicitly 
fitted to do so. On the other hand, the predictions pro- 
vided by this model do not reach the quantitative accu- 
racy of the results obtained in the fixed-cell case. More 
precisely, Table IIVI shows significant differences between 
the first-principles results (labeled "LDA" ) and the pre- 
dictions of the L° model for the structure of the tetrag- 
onal ground state, especially as regards the aspect ratio 
(c/a) of the unit cell and the participation of the Pb 
atoms in the ferroelectric distortion. 

Wanting to increase the model's accuracy, we decided 
to improve the description of the strain-phonon couplings 
by adding the SAT represented by (Pb x — 02 x ) 2 r] 1 , where 
we use the compact notation introduced above (see Fig. [3] 
and Table ITT|) . Note that the resulting model, which we 
label L 1 , combines A-like terms, whose values are fixed to 
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TABLE IV. Results for the tetragonal ground-state structure of PTO with relaxed cell parameters. Lattice vectors and atomic 
displacements are given in Angstroms. The atomic displacements are as described in the caption of Table IIIII We show the 
first-principles results (first row) followed by the results obtained from models with different descriptions of the strain-phonon 
coupling terms (see the text). Energies are given in meV/f.u., and we take Ens as the zero of energy. 



Method 


a 


c 


c/a 


UPbz 


UTiz 


UOlz 


U03z 


Energy 


LDA 


3.864 


3.974 


1.029 


0.230 


0.106 


-0.133 


-0.071 


-37.7 


model L° 


3.908 


3.987 


1.020 


0.200 


0.103 


-0.122 


-0.060 


-34.5 


model L 1 


3.863 


3.968 


1.027 


0.220 


0.099 


-0.128 


-0.063 


-39.9 


model L" 


3.861 


3.978 


1.030 


0.227 


0.102 


-0.132 


-0.066 


-43.1 


model L UI 


3.856 


3.968 


1.029 


0.221 


0.098 


-0.128 


-0.062 


-39.9 



(a) model L° (b) model L 1 (c) model L 
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FIG. 10. (Color on line.) Force-constant bands corresponding to PTO structures with m = and subject to an uniaxial 
strain of 2% (stretching along z), as obtained from three models that include different strain-displacement couplings (see the 
text). Black solid lines are the model results, and red dashed lines depict the bands obtained from first principles. Note the 
non-analytic behavior of the spectrum when approaching T from the [(fa, 0,0] or [0, 0, q z ] directions. 



those of the L° potential, with one A-like free adjustable 
parameter. Such a parameter was fitted to better repro- 
duce the ground state structure; as shown in Table IIV1 
this led to a significant improvement over the L° result. 
Further improvement of the c/a value can be achieved 
by additionally introducing the higher-order SAT repre- 
sented by (Pb x — 02 x ) 2 rj{ (model L 11 ), at the expense of 
worsening the agreement for other structural parameters 
and energies. 

Let us mention here another model-construction ex- 
periment that we made. Noting the importance of the 
strain-phonon couplings in PTO, one may wonder which 
are the interaction terms responsible for the main ef- 

( 1 2") 

fects. By inspecting the ' ; parameters computed 
directly from first principles, it is easy to identify the 
two most prominent ones, which involve Pb-0 and Ti-0 
nearest-neighboring pairs. More specifically, the key cou- 
plings are captured by the A-like parameters number 2 
and number 13 from Table |TTJ Hence, we considered a 
model that includes only these two strain-phonon cou- 
plings (L 111 ); interestingly, as shown in Table llVl such 
a simple potential is able to render good results for the 
structure and energy of PTO's ground state. 

The quality of these models can be further evaluated 
by checking how well reproduce the first-principles results 
for the force-constant bands of strained configurations. 
As already mentioned, Fig. [TO] shows an essentially per- 
fect agreement for model L° , which is largely preserved in 
models L 1 and L 11 (the latter is not shown). Naturally, 



the agreement is worse for the minimal model L . Fig- 
ure [TT] also shows the results that model L 1 gives for the 
force-constant bands of PTO's tetragonal ground state, 
as compared with the first-principles calculations. As 
in the fixed-cell cases of Fig. [8l it is apparent that the 
considered model is not sufficient to render a precise de- 
scription of all the bands. Yet, the qualitative agreement 
is satisfactory. 



4- Temperature- dependent behavior 

We studied the T-dependent behavior of our PTO 
models by running MC simulations in which both the 
atomic displacements and strains were allowed to ther- 
mally fluctuate. Figure [T2l shows the basic results for our 
L° model when simulated in two different situations: (i) 
under the condition of zero external pressure and (ii) by 
imposing an external hydrostatic pressure of —13.9 GPa, 
which counteracts the underestimation of the LDA re- 
sult for the cubic lattice constant. (Taking as a reference 
the cubic lattice constant obtained by extrapolating to 
K the experimental results in Ref. [42], this underesti- 
mation can be approximated to be about 2.2%.) Note 
that this kind of correction is customarily made in LDA- 
bascd effectivc-Hamiltonian works,— ~— and we adopt it 
here for the sake of an easier comparison with the litera- 
ture. 

As can be appreciated in Fig. [T^J our simulated PTO 
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FIG. 11. (Color on line.) Force-constant bands for the ground 
state of PTO. Black solid lines depict the results from model 
L 1 , and red dashed lines show the results obtained from first 
principles. 

undergoes a phase transition from the high-T cubic phase 
to a low-T tetragonal structure in which one polarization 
component (z in our default Cartesian setting) becomes 
different from zero. The transition is accompanied by a 
deformation of the cell, which acquires a c/a > 1 aspect 
ratio. The computed Curie temperature is about 225 K 
when no external pressure is applied, and increases to 
about 450 K when we correct for the LDA overbinding. 
This is the expected behavior, as it is known that the 
strength of the FE instabilities in these perovskite oxides 
is very sensitive to volume changes (which is the reason 
why they have very good piezoelectric properties). 

Figure [T3] shows the results obtained for all the model 
potentials listed in Table IIVI simulated under the same 
hydrostatic pressure of —13.9 GPa. Remarkably, in spite 
of their similarly good description of the ground state en- 
ergy and structure, we observe very large differences in 
the predicted Tc's. It is interesting to note that, contrary 
to what we would have expected ) 52 ' 53 the obtained Tc's 
do not correlate well with the energy difference between 
the ground state and the RS, nor with the magnitude 
of the FE distortion. Thus, for example, the lowest Tc 
(about 440 K) corresponds to the L 111 potential, in spite 
of the fact that the weakest FE instability (c/a = 1.020; 
E gs — Ens = —34.5 mcV/f.u., where E gs is the ground 
state energy) corresponds to the L° model. (The same 
trends were observed in the MC runs with no applied 
pressure.) It is thus clear from these results that the 
computed Tc's are strongly dependent on details of the 
PES that are not reflected in the energy and structure of 
the ground state, a conclusion that can be extended to 
all physical properties that we may obtain from our MC 
simulations. Hence, the results in Fig. [13] evidence the 
critical importance of developing models that include all 
the atomic degrees of freedom, and allow for a system- 
atic improvement of the PES description, if we want to 
obtain accurate first-principles results of the thermody- 
namic properties of materials like PTO. 

Let us conclude by giving an additional and striking 
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FIG. 12. (Color on line.) Temperature-dependent polariza- 
tion [panel (a)] and strains [panel (b)] of PTO as calculated 
from our L° potential in two different conditions of external 
pressure (see the text). The LDA-relaxed cubic structure de- 
fines the zero of strain. The value of the thermal expansion 
coefficient (a) of the high-temperature phase is indicated. 



example of the importance of hidden atomistic effects in 
determining the macroscopic properties of this material. 
Our best model for PTO is probably the one labeled L 1 , 
which renders a FE transition at Tc « 510 K. Interest- 
ingly, Waghmare and Rabe (WR) constructed an effec- 
tive Hamiltonian for PTO, considering only polar local 
modes and strains as the model variables, that results 
in a significantly higher Tc of about 660 Kp^i At first 
sight such a discrepancy may seem surprising, and we 
made an effort to understand its origin in some detail. 
First, we checked that our model reproduces the ener- 
getics of the FE instabilities given by the WR Hamil- 
tonian rather closely, despite the differences in the first- 
principles calculations (e.g., in the pseudopotentials) em- 
ployed to compute the parameters. Further, we ran simu- 
lations with modified versions of our model to test subtle 
features of the WR energy parametrization (e.g., the in- 
clusion of high-order terms for the polar local modes), 
and concluded that they cannot account for the discrep- 
ancy in the computed Tc. 

We thus turned our attention to the qualitatively dis- 
tinct features of our model. Most notably, we describe 
not only the FE instabilities and strains, but also the un- 
stable AFD distortions sketched in Fig. [5] It is known 
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FIG. 13. (Color on line.) Same as Fig. 1121 but for the refined 
PTO model potentials discussed in the text. In all cases an 
external pressure of —13.9 GPa is applied. In panel (a) dashed 
lines show the result for L°. 



that, in most perovskite oxides, the interaction between 
FE and AFD modes is a competitive one, so that they 
tend to suppress each other.— Hence, to evaluate the ef- 
fect of such a competition in our simulated PTO, we ran 
simulations in which the Oq rotational modes were not al- 
lowed. We imposed this constraint by restricting the mo- 
tion of the oxygen atoms as shown in the sketch of Fig. [Til 
Let us stress that such a constraint docs not affect the 
energetics associated with the development of the spon- 
taneous polarization, the FE ground state being exactly 
retained. Figure shows the results for our L 1 model: 
In the case without AFDs we got Tq w 825 K, which 
lies about 300 K above the result obtained from the un- 
constrained simulation. (The FE-AFD competition was 
predicted by other authors to have a similarly large im- 
pact on the Tq of the PhZri^Ti^Oa solid solution^) 
The details of these FE-AFD interactions will be dis- 
cussed at length in a future publication. We show this 
result here just as a striking example of the physical ef- 
fects that we are likely to miss if we restrict ourselves to 
effective models that, in spite of looking complete (as e.g. 
they capture the basic features of the ground state), may 
turn out to be too simple. 
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FIG. 14. (Color on line.) Temperature-dependent polariza- 
tion as obtained from MC simulations of our L 1 model. The 
black solid symbols show the results obtained when we allow 
all possible atomic movements (as pictorially depicted for the 
oxygen atoms in the left inset); the red open symbols show 
the results when we suppress the oxygen displacements asso- 
ciated with the rotations of the Og octahedra (the right inset 
shows the allowed oxygen displacements in this case). 



C. SrTiOs 

1. Harmonic terms Eha,r({ui}) and E s (rj) 

We extracted all the non-zero harmonic coupling terms 
from DFPT calculation s 30 ! 35 carried out with Abinit4°- 
Representative results for the force-constant bands of the 
cubic RS are shown in Fig. [Jj)] Note that in this case 
we only have AFD-related instabilities; indeed, our LDA 
calculations render low-energy, but perfectly stable, FE 
modes for the cubic phase of STO. 

The short-range interatomic interactions that we ob- 
tained for STO have the same spatial extent as those 
computed for PTO and described above. As regards the 
electronic dielectric tensor, Born effective charges, and 
harmonic elastic constants, our results for STO are also 
analogous to the ones for PTO described above, as the 
number of symmetry-independent terms is the same for 
both materials. 



2. Fitting E^ nh ({ui}) 

We fitted the terms in £7 a nh ({ u i}) by working with a 
training set of relevant low-symmetry phases that main- 
tain the cubic STO cell (i.e., with rj = 0). More pre- 
cisely, we considered the following AFD-distorted struc- 
tures: AFD* (PA/mbm space group), AFD" (7 4/ 'mem), 
AFD" Z (Imma, with rotations of equal amplitude about 
x and z) and AFD° yz (i?3c). As in the case of PTO, we 
determined such low-symmetry structures ab initio by 
distorting the RS according to a specific unstable eigen- 
vector and relaxing the resulting structure while preserv- 
ing the targeted symmetry. The energies and distortion 
amplitudes computed for these structures are given in 
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FIG. 15. (Color on line.) Same as Fig.fl] but for SrTi0 3 . 

Table El 

Additionally, in our set we also included two structures 
generated by distorting the cubic phase according to tin- 
lowest-energy FE eigcndisplacement obtained from our 
DFPT calculations [which strongly resembles the typi- 
cal FE unstable mode depicted in Fig. [SJa)]. More pre- 
cisely, we considered two distortions involving polariza- 
tions along the [001] and [111] Cartesian directions, re- 
spectively. We included such structures in the training 
set to capture the anharmonicity of the low-lying FE 
modes, which should play a role in determining the non- 
linear dielectric response properties of interest in STO. 

As in the case of PTO, we worked with a relatively 
simple model restricted to pairwise anharmonic interac- 
tions extending up to the 4th order of the Taylor series. 
Further, we restricted ourselves to interactions between 
first-nearest-neighboring atoms, an approximation that 
is justified by the rapid spatial decay of the anharmonic 
corrections that we observed for STO as well. As we 
have mentioned already, these truncations result in the 
15 SATs listed in Table H 

Additionally, in the case of STO we tried to identify 
the minimal set of SATs that capture the energetics of the 
low-symmetry structures in our training set. We found it 
possible to do so by considering only the 4th-order terms 
with numbers 4, 12, 13, and 15 in Table Q] We computed 
the corresponding parameters by optimizing QT e and 
Q^VEi and obtained the results summarized in Fig. 1161 
and Table IVl The agreement with the first-principles data 
is very good, and we checked that no significant improve- 
ment is obtained by including other SATs listed in Ta- 
ble n 

To construct the present model of STO, we did not 
take the additional step of minimizing a goal function 
GJ°hoss with information about the Hessian matrices of 
the low-energy structures. It is therefore interesting 
to check whether the mode stiffnesses calculated using 
our effective potential reproduce well the first-principles 
data. Representative results are depicted in Fig. IPTTa). 
where we show DOS plots constructed from the force- 
constant eigenvalues X qs for the AFD" Z structure, which 
is the predicted ground state of the material for rj = 
(see Table |V| . As it can be seen, our model properly de- 
scribes the structure as being a stable one (i.e., we find 
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FIG. 16. (Color on line.) Potential-energy wells connecting 
the RS of STO with the low-symmetry phases mentioned in 
the text. The results obtained from our model potential are 
shown with lines, and the points indicate the first-principles 
results for the energy minima or saddles, as well as for the 
two FE-distorted structures considered. All the states shown 
preserve the cubic cell of the RS (rj — 0). The AFD" Z and 
AFD" yj; curves are essentially on top of each other and cannot 
be distinguished. 



no modes with X qs < 0) and reproduces well the general 
shape of the spectrum. In fact, the overall agreement for 
the spectrum of Hessian eigenmodes is comparable to the 
one obtained for PTO. 



3. Strain-phonon term E ap ({ui},r}) 

As in the case of PTO, we considered only the lowest- 
order terms that are allowed by symmetry and cap- 
ture the most important strain-phonon effects, which are 
given by the coefficients A*- 1 ' 2 -*. We computed them di- 
rectly by employing the finite-difference approach sum- 



(a) AFD^ (77 = 0) 



(b) AFDJ (relaxed) 




DOS (LDA) DOS (model) 



DOS (LDA) DOS (model) 



FIG. 17. Density of states (DOS) for two representative low- 
energy phases of STO, obtained in a way that is analogous 
to the one described in the caption of Fig. [4] We show the 
results obtained from first principles and from the presented 
model potential. 
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TABLE V. Structural parameters of STO's low-energy phases (see the text) as obtained from first-principles LDA calculations 
and from the presented model. A cubic cell with a = 3.845 A was used in all cases, except in the one marked with an asterisk; in 
that case, a full structural relaxation was performed and the resulting pseudocubic lattice constants are given. The amplitude 
of the AFD modes is quantified by the corresponding C>6 rotation angle given in degrees (in the AFD!J Z and AFD" HZ cases, 
we have equal-magnitude rotations about two and three Cartesian axes, respectively; we give the rotation angle for one axis.). 
Energies are given in meV/f.u., taking the result for the RS as the zero of energy. Note that the AFD" Z structure displays 
additional small distortions; for example, there are anti-polar displacements of the Sr atoms, the off-centering being about 
0.008 A and about 0.002 A for the LDA and model calculations, respectively. 



Structure 



Method Energy 06 rot. 



AFD" mem) 
AFD° Z (Imma) 



AFD" yz (R3c) 



AFD Z (PA/mbm) 
AFD" (7 4 I mem) * 



LDA 


-18.9 


6.7 


model 


-18.8 


6.6 


LDA 


-19.4 


4.9 


model 


-20.0 


4.9 


LDA 


-18.8 


3.9 


model 


-19.7 


3.9 


LDA 


-3.7 


4.5 


model 


-3.2 


4.3 


LDA 


-23.0 


7.5 


model 


-23.0 


7.4 



3.825 A 
3.824 A 



3.869 A 
3.867 A 



marized by Eqs. (|32|) and (|38|) . using strains of ±2% for 
the finite-difference calculations. The resulting model de- 
scribes the correct ground state of STO, which is char- 
acterized by an AFD" distortion (14/ mem space group). 
Note that, for r\ = 0, our first-principles calculations in- 
dicate that the ground state is associated with an AFD° Z 
distortion. Hence, as in PTO's case, the strain-phonon 
couplings play a key role in determining the symmetry 
of the lowest-energy structure; also like in the case of 
PTO, such an effect is captured by the A^ 1 ' 2 -* parameters 
computed directly via our finite-differences scheme, even 
though they were not explicitly fitted to do so. 

This model gives an excellent quantitative description 
of STO's ground state (see Table IV]) . and reproduces 
reasonably well the corresponding Hessian matrix [see 
Fig- ITWb)]. Hence, we took this potential as our effec- 
tive model for STO, without any further refinement. 



4- Temperature- dependent behavior 

Figure [18] shows the basic results from the MC simu- 
lations performed with our model potential for STO. As 
in the case of PTO, we ran simulations (i) under zero 
applied pressure and (ii) under an expansive hydrostatic 
pressure of —9.2 GPa, which approximately corrects for 
the LDA overbinding. [To compute the correction, we 
used as reference a cubic lattice constant of 3.90 A, ob- 
tained by extrapolating to K the experimental results 
for the cubic phase of STO in Ref. In both cases we 
get a phase transition from the high-T cubic phase to a 
low-T tetragonal structure (1 4 /mem) with AFD" charac- 
ter. The transition temperature is about 300 K when no 
pressure is applied, and decreases to about 160 K upon 
application of —9.2 GPa. This is the expected behav- 
ior, as the applied pressure is known to (i) reduce the 
strength of the AFD instabilities and (ii) enhance the 



FE-AFD competition by softening the FE distortions. 
These effects have been studied in previous theoretical 
works on STO and related perovskitesj 13 i 55 we have cap- 
tured them implicitly (i.e., without any ad hoc fitting) 
when constructing our model. 

D. Discussion 

Let us conclude this Section by commenting on how 
well our models reproduce experiment; in particular, let 
us focus on their performance to predict one of the most 
basic properties of these ferroic materials, namely, the 
temperature of their structural transition. 

The cffcctivc-Hamiltonian approach to FE pcrovskitcs 
has been very successful in reproducing non-trivial be- 
haviors of many complex materials qualitatively; exam- 
ples include the phase diagram of chemically-disordered 
solid solutions ; 15 ^ 6 the occurrence of multiferroic 
orders ^ 18 strain—, finite-size— ~—, and electrostatics- 
58 i 59 driven effects, and quantum-phase transitions! 14 i 60 
However, whenever the model parameters have been ob- 
tained directly from first principles, and despite the use 
of pressure corrections as the one considered here, the 
quantitative agreement for the predicted transition tem- 
perature has been a poor one. Thus, for example, the 
cubic-to-tetragonal transition of BaTiOs was predicted 
to occur at about 300 while the experimental re- 

sult is 400 K. In the case of KNbOa^i the simulations 
render a cubic-to-tctragonal transition at 370 K, while 
the experimental Tq is about 700 K. In this context, the 
result of Waghmare and Rabeii for PTO - i.e., a Tq of 
660 K that compares reasonably well with the observed 
value of 760 K - might be considered as an example of 
good agreement between theory and experiment. 

The difficulties of the cffcctivc-Hamiltonian method to 
obtain correct transition temperatures were analyzed in 
Ref. |H, where it was suggested that the discrepancy is to 
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FIG. 18. (Color on line.) Temperature-dependent AFD a dis- 
tortions [panel (a)] and strains [panel (b)] of STO as obtained 
from MC simulations of our model. The LDA-relaxed cubic 
structure defines the zero of strain. The value of the ther- 
mal expansion coefficient (a) of the high-temperature phase 
is indicated in panel (b). 



be partly attributed to an incorrect description of ther- 
mal expansion, a problem that is a direct consequence 
of the coarse-graining step involved in the construction 
of the model. Indeed, the effective Hamiltonians tend to 
give an essentially null thermal ex pan sion at high tem- 
peratures (see e.g. Fig. 10 in Ref. [ill ), which is clearly 
against the experimental evidence. 

Our models take into account all the degrees of free- 
dom in the material, and we can thus hope to improve 
on this aspect. As can be seen in Fig. Q21 for PTO 
we get a thermal expansion coefficient at high temper- 
atures between 8.3xl0~ 6 K _1 (from simulations with 
no applied pressure) and 9.1 xlO -6 K" 1 (obtained when 
pressure is applied to correct for the LDA overbinding) , 
to be compared with the experimental result of 12.6 x 
10~ 6 K -1 42 Similarly, Fig. [TBI shows a thermal expan- 
sion between 7.2xlCT 6 K" 1 and 8.3xl0" 6 K" 1 for STO 
at high temperatures, to be compared with the value 
of 8.8xl0 -6 K _1 obtained from experiments^ Hence, 
our models clearly improve the effective-Hamiltonian de- 
scription of this effect; yet, the discrepancy between the 
computed and measured transition temperatures remains 
present. 



In the case of PTO, we have found solid evidence that 
Tq depends very significantly on details of the PES that 
are most often ignored in theoretical works. In particu- 
lar, our results strongly suggest that a realistic model for 
PTO must necessarily include both FE and AFD degrees 
of freedom, as their competition is far from being negligi- 
ble. Accordingly, we should probably consider as partly 
fortuitous the relatively accurate result obtained for Tq 
in Ref. [TJ where a model without AFDs was employed. 

Our results lead to the following important conclusion: 
We cannot expect to obtain accurate values for PTO's Tc 
from models that only reproduce the basic first-principles 
results (i.e., energy and structure) for the low-symmetry 
phases of the material. Further, in order to improve 
the agreement with experiment, we should probably ex- 
tend our model to better reproduce the lattice-dynamical 
properties of the key low-energy structures and other de- 
tails of the PES. Our results in Fig. [TBI suggest that im- 
provements of that sort, even though they may look like 
second-order corrections to the relevant PES, can actu- 
ally affect the computed Tq by as much as 100 K. 

On the other hand, such a strong sensitivity to the 
details of the potential has important implications re- 
garding the accuracy required from the first-principles 
methods used to compute the model parameters. To- 
gether with the incorrect treatment of thermal expansion, 
the authors of Ref. mentioned DFT inaccuracy as the 
second reason to explain the fact that the Curie temper- 
atures obtained from effective-Hamiltonian simulations 
are typically too low as compared with the experimen- 
tal ones. Their conjecture was that the FE instabilities 
obtained from first principles were too weak, meaning 
that DFT was probably underestimating the \E gs — -Ers| 
energy difference between the RS and the FE ground 
state. Our results show that, while important, this is 
by no means the only characteristic of the PES that has 
a large impact on the computed To Hence, to get ac- 
curate results, we need a first-principles theory that not 
only describes correctly the energetics of the FE insta- 
bility, but also captures accurately more subtle PES fea- 
tures such as the anharmonic couplings between different 
structural distortions, including those that do not partic- 
ipate directly in the transitions. These are very demand- 
ing requirements for our simulation techniques; thus, it is 
unclear whether we presently count with first-principles 
methods that can predict an accurate Tc for materials 
like PTO. 

Most of the above considerations probably apply to 
STO as well. Yet, as the transition occurs at a rela- 
tively low temperature in this case, an additional fac- 
tor must be taken into account. As demonstrated in 
a variety of theoretical works ; 14 ' 60 ' 63 in order to get a 
precise calculation of the structural transition tempera- 
tures of ABO 3 perovskites, it is important to consider 
the quantum (i.e., wave-like) character of the atoms. In- 
deed, Zhong and Vanderbilt simulated STO at both the 
classical and quantum-mechanical levels, using an effec- 
tive Hamiltonian constructed from first principles, and 
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found that quantum fluctuations shift down the cubic- 
to-tetragonal transition by about 20 KJ^ (Quantum ef- 
fects will typically promote disorder, and thus result in 
reduced transition temperatures.) Hence, in the case of 
STO, we can assume that part of the discrepancy be- 
tween our computed transition temperature (~ 160 K) 
and the experimental one (105 K) comes from the fact 
we treated atoms as classical objects in our MC simula- 
tions. Finally, let us note that our computed transition 
temperature seems consistent with the result of about 
130 K reported in Ref. [ID for the classical case. 

IV. SUMMARY AND CONCLUSIONS 

We have proposed a new method for the construction 
of first-principles model potentials that permit large-scale 
simulations of lattice-dynamical phenomena. Our scheme 
mimics the traditional approach to lattice dynamics in 
solid-state textbooks, i.e., we start from a suitably cho- 
sen reference structure (RS) and express the energy of the 
material as a Taylor series for the structural distortions 
of such a RS. There are many advantages in adopting 
such a simple approach; most notably, our potentials can 
be trivially formulated for any compound, and their abil- 
ity to reproduce the first-principles data can be improved 
in a systematic and well-defined way. Further, most of 
the potential parameters correspond to the usual (elas- 
tic, force-constant, etc.) tensors discussed in condensed- 
matter theory, which allows for a transparent physical 
interpretation. 

We have described the details of such an approach, 
and proposed a practical strategy to compute the model 
parameters from first principles. Our method is espe- 
cially convenient in that regard too, as we have shown 
that many of the key model parameters (e.g., all the cou- 
plings at the harmonic level) can be readily obtained from 
density-functional-pcrturbation-theory calculations that 
are widely available today. 

We have illustrated our method with applications to 
two especially challenging cases, namely, ferroic oxides 
PbTiOs and SrTi03. These materials undergo struc- 
tural phase transitions driven by soft phonon modes, 
which implies that the potential-energy surface (PES) 



that our models have to capture is strongly anharmonic. 
We have discussed in detail the case of PTO, where 
the large structural deformations involved in the ferro- 
electric phase transition make it especially challenging 
to construct a quantitatively accurate model. More- 
over, we have solved our PTO potential by means of 
Monte Carlo simulations and discovered a variety of un- 
expected effects, ranging from novel structural phases 
when the strain deformations are constrained to a sur- 
prisingly strong dependence of the computed Curie tem- 
perature on the details of the PES. The case of STO 
turned out to be much easier to tackle and led to quanti- 
tatively more accurate predictions, probably because the 
structural distortions involved in its ferroic transforma- 
tion are smaller. The connections of our method with the 
so-called Erst-principlcs effective Hamiltonian approach 
to the study of temperature-driven effects in ferroic per- 
ovskitc oxides - which was introduced about 20 years 
ago, and of which our scheme can be considered a nat- 
ural extension and generalization - have been discussed 
in some detail. 

We believe that our effective potentials can be used 
to great advantage in the investigation of the thermody- 
namic properties of (meta)stable material phases, which 
are largely dominated by harmonic effects that our mod- 
els describe with first-principles accuracy. While we have 
not considered any such case in this work, we believe 
that the demonstrated ability of our models to deal with 
strongly anharmonic effects suggests that their applica- 
tion to (the much simpler) quasi-harmonic cases will be a 
very successful one. Hence, we hope the current method- 
ology can become a standard tool for large-scale simula- 
tions of the lattice-dynamical properties of materials at 
realistic operating conditions of temperature, pressure, 
etc. 
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